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□  □ 


Abstract 


The  goal  of  this  dissertation  is  to  develop  a  non  parametric  method  for  obtain¬ 
ing  a  confidence  interval  for  the  mean  of  a  stationary  sequence.  As  indicated  in 
the  Literature,  nonparametric  confidence  intervads  in  practice  often  have  undesir¬ 
able  small-sample  asymmetry  and  coverage  characteristics.  These  phenomena  are 
partially  due  to  the  fact  that  the  third  and  fourth  cumulants  of  the  point  estima¬ 
tor  for  the  stationary  mean,  unlike  those  of  the  standard  normal  random  variable, 
are  not  zero.  We  will  apply  Edgeworth  and  Cornish- Fisher  expansions  to  obtain 
asymptotic  expansions  for  the  errors  associated  with  confidence  intervals.  The 

analysis  isolates  various  elements  that  contribute  to  errors  and  makes  it  possible 

0 

for  us  to  estimate  each  element  and  hopefully  correct  the  errors  to  a  smaller  order. 
We  will  use  Glynn's  method  to  develop  first  and  second  order  pivots  for  the  confi¬ 
dence  intervals.  Furthermore,  these  procedures  also  improve  the  asymptotic  order 
of  confidence  interval  accuracy.  7  ^  ' 

Keywords:  simulation,  steady-state,  nonparametric  confidence  intervals, 
Edgeworth  expansions,  stationary  processes 
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Chapter  1 


Introduction 


The  use  of  stochastic  analysis  has  long  been  an  important  tool  in  the  study  of 
complex  systems  such  as  manufacturing,  computer,  and  communications  systems. 
As  the  complexity  of  systems  grow,  it  appears  to  be  the  rule  rather  than  the  excep¬ 
tion  that  many  detailed  stochastic  models  are  so  complex  such  that  it  is  extremely 
difficult  or  impossible  to  obtain  an  exact  analytic  solution.  Although  simulation 
was  often  viewed  as  a  “method  of  last  resort”  to  be  used  only  when  everything  else 
failed,  the  intrinsic  complexity  of  stochastic  models  for  real  world  systems,  recent 
advances  in  simulation  methodology,  and  the  availability  of  computing  power  and 
software  packages  have  made  simulation  one  of  the  most  widely  used  tools  in  system 
analysis  and  operations  research. 

In  a  typical  simulation  study,  we  first  construct  a  mathematical  model  that 
corresponds  to  the  real  world  system  to  be  studied.  We  then  perform  computer 
sampling  experiments  on  the  model,  collect  and  analyze  the  output  sequences,  and 
make  inferences  about  the  behavior  of  the  system.  Since  simulation  is  a  statistics 
based  computer  sampling  procedure,  appropriate  statistical  techniques  must  be 
used  in  both  design  and  analysis  of  the  simulation  study  before  any  meaningful 
conclusion  can  be  drawn. 

Simulation  output  analysis  is  an  area  of  active  research.  In  most  situations, 
we  want  to  estimate  quantities  associated  with  the  stochastic  system  being  simu¬ 
lated.  When  estimating  a  quantity  only  by  a  sample  mean,  no  indication  of  the 
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variability  of  the  estimate  is  given.  The -standard  deviation  of  the  estimate  is  one 
device  used  for  indicating  the  reliability  or  precision  of  an  estimate.  What  is  more 
informative,  however,  is  a  confidence  interval.  Basically,  we  first  obtain  a  point  es¬ 
timate,  calculate  its  standard  deviation,  and  then  construct  a  confidence  intervaJs 
for  each  quantity  to  access  the  statistical  variability  of  the  point  estimate.  It  is 
generally  thought  that  the  construction  of  confidence  interv£ds  is  one  of  the  most 
difficult  and  important  issues  in  the  field  of  simulation  output  analysis. 

Although  there  is  a  large  statistical  literature  devoted  to  the  assignment  of 
confidence  intervals,  most  of  it  can  not  be  directly  applied  to  simulation  output 
analysis.  One  important  reason  for  this  is  that  most  stochastic  processes  associ¬ 
ated  with  real  world  sim  ilation  studies  do  not  satisfy  the  standard  assumptions  in 
the  statistical  literature,  namely,  independence,  stationarity,  etc.  At  present,  there 
are  many  methods  proposed  in  the  simulation  literature.  The  most  common  meth¬ 
ods  are  replications,  batch  means,  overlapping  batch  means,  regenerative,  spectral, 
autoregressive,  autoregressive  moving  average,  and  standardized  time  series. 

On  the  other  hand,  the  recognition  that  computing  costs  are  decreaising  has 
allowed  statisticians  and  simulation  researchers  to  consider  confidence  intervals 
methods  which  are  computationally  more  intensive  but  statistically  better  behaved 
than  previous  techniques.  Among  these  methods  are  the  jackknife,  bootstrap,  other 
resampling  plans,  and  Johnson-Glynn  pivital  tranformations. 

The  underlying  mathematical  model  in  this  study  is  a  stationary  stochastic 
process  which  satisfies  certain  regularity  conditions.  As  in  any  steady  state  analysis 
of  time  series,  a  key  assumption  is  that  the  initial  state  should  have  a  very  small 
effect  on  the  overall  behavior  of  the  system.  However,  when  studying  a  specific 
real  world  system,  this  assumption  does  not  always  hold.  On  the  other  hand. 
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initial  effects  decay  exponentially  over  time  for  several  classes  of  stochastic  processes 
so  that  steady-state  simulation  is  often  appropriate  for  long  simulation  runs.  In 
Chapter  5,  we  give  several  numerical  examples  showing  that  their  initial  conditions 
do  not  affect  the  behavior  of  the  confidence  intervals. 

In  this  dissertation  we  develop  a  nonparametric  method  for  obtaining  a  more 
accurate  asymptotic  confidence  intervaJ  for  the  sample  mean  of  a  stationciry  process 
which  satisfies  certain  regularity  conditions.  By  “more  accurate  asymptotic  confi¬ 
dence  interval”  we  mean  that  the  coverage  error,  which  is  the  difference  between 
the  nominal  and  the  actual  coverage,  associated  with  our  confidence  interval  will 
be  of  lower  order  than  that  of  the  traditional  methods. 

Oar  starting  point  is  the  traditional  batch  mecins  method.  We  use  the  idea 
of  Johnson-Glynn  pivotal  transformations  to  obtain  better  confidence  intervals  for 
the  quantities  of  interest.  No  assumption  has  been  made  that  the  observed  data  are 
sampled  from  either  i.i.d.,  regenerative,  or  ARMA  processes.  We  believe  this  more 
general  case  is  robust  for  the  output  analysis  of  real  world  simulation  experiments. 
The  procedures  we  propose  do  not  require  the  selection  of  any  critical  constants 
that  can  not  be  reasonably  preset.  This  fact  and  the  less  restrictive  nature  of  the 
assumptions  will  also  allow  us  to  implement  these  procedures  as  a  software  package 
for  the  output  analysis  of  many  real  world  simulation  studies. 

The  baisic  approach  of  the  procedures  we  propose  follow  from  Johnson  [15]. 
Glynn  [12],  and  Titus  [25].  As  indicated  in  Glynn  [12],  nonparametric  confidence 
intervals  in  practice  often  have  undesirable  small-sample  asymmetry  and  coverage 
characteristics.  We  will  apply  Edgeworth  expansion  theory  to  obtain  asymptotic 
expansions  for  the  errors  associated  with  confidence  intervals.  The  analysis  isolates 
the  various  elements  that  contribute  to  the  errors.  We  then  use  Glynn's  method  to 
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develop  first  aiid  second  order  pivots  to  the  confidence  intervaJs,  which  deal  with 
asymmetry  problems  and  coverage  difficulties,  respectively.  These  procedures  also 
improve  the  asymptotic  order  of  confidence  interval  accuracy  in  the  sense  that  the 
actual  coverage  of  the  corrected  confidence  interval  is  closer  to  the  nominal  coverage 
rate  than  previous  methods. 

Johnson  [15]  is  perhaps  the  first  author  to  use  these  procedures.  He  derives 
a  first  order  pivot  for  the  t-statistic  for  independent  and  identically  distributed 
samples.  Glynn  [12]  extends  this  idea  to  a  second  order  pivot  for  the  ratio  estimators 
of  the  regenerative  processes.  Titus  [25]  applies  the  same  idea  to  cisymptotically 
stationary  autoregressive  processes  of  finite  order. 

The  organization  of  this  dissertation  is  as  follows.  In  Chapter  2,  we  develop  the 
necessary  background  for  cumulants,  Edgeworth  expansions,  and  Cornish-Fisher 
expansions.  Some  properties  of  the  cumulants  for  a  stationary  process  satisfying 
some  regularity  conditions  are  also  discussed  there.  In  Chapter  3,  we  derive  the 
first  and  second  order  Johnson-Glynn  pivotal  transformations  to  correct  the  error 
in  confidence  interval  coverage.  Some  computational  and  theoretical  issues  are 
also  discussed  there.  To  demonstrate  how  our  method  works,  several  numerical 
examples  are  displayed  in  Chapter  4.  Chapter  5  summ2irizes  the  strength  and 
weaJcness  of  the  Johnson-Glynn  pivots  as  applied  to  batch  means  method.  The 
Appendix  contains  some  of  the  more  technical  proofs. 


Chapter  2 


Mathematical  Background 


In  this  chapter,  we  will  review  the  necessary  mathematical  background  re¬ 
quired  for  this  dissertation.  In  Section  2.1,  definition  of  cumulants  as  well  as  some 
important  properties  of  cumulants  will  be  given.  Section  2.2  concerns  the  Edge- 
worth  expansion  while  the  Cornish- Fisher  expansion  is  discussed  in  Section  2.3.  In 
Section  2.4,  some  new  uniqueness  properties  of  the  Cornish-Fisher  expansion  will 
be  presented.  Section  2.5  studies  some  properties  associated  with  the  cumulants 
of  stationary  processes.  Finally,  in  Section  2.6,  we  will  present  results  about  the 
consistency  properties  of  the  sample  moments  for  stationary  processes. 

2.1.  Cumulants 

Consider  a  reindom  variable  X  with  distribution  function  .  The  moments 
fi'r  =  EX^,  r  >  0,  and  central  moments  pr  =  E{X  —  EXy,  r  >  0,  are  useful 
constants  for  measuring  properties  of  X  and,  in  some  cases,  uniquely  characterize 
the  distribution  function  (Chung  [8],  pp.  98-99).  For  some  statistical  analyses 
another  sequence  of  constcints,  the  cumulants,  are  more  useful  from  a  theoretical 
point-of-view. 

Let  be  the  characteristic  function  of  X.  Then,  subject  to  existence,  the 
cumulants  of  X,  k/s,  are  formally  defined  by  (Kendall  and  Stuart  [16],  p.  69) 

^xll)=  r 

J — oo 
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r=0 


r=l 


(2.1) 


Thus  Kr  is  the  coefficient  of  [Uy jr\  in  log$;f(-)  when  the  expansion  in  power  series 
exists;  log$x  is  called  the  cumulant  generating  function  (c.g.f.)  or  the  second 
characteristic  function.  However,  the  first  terminology  is  somehow  misleading  in 
the  sense  that  log  exists  in  some  neighborhood  of  the  origin  even  if  the  moments 
and  cumulants  do  not  exist. 


The  proceeding  expansions  can  be  made  rigorous.  It  is  known  that  $  v  is 
uniformly  continuous  in  (—00,00)  and  4>x(0)  =  1  (Chung  [8],  p.  143).  Therefore 
there  exists  a  neighborhood  of  the  origin  in  which  is  different  from  zero;  let 
(f(  <  A  be  this  neighborhood.  By  taking  the  principal  branch  of  the  logarithm. 
log^>;f  can  be  uniquely  defined  for  |t|  <  A.  Moreover,  log<^;f  is  continuous  and 
vanishes  at  t  =  0.  Now  assume  for  some  r  >  1,  E\XY  exists,  we  then  have 

=  (2-2) 
j=Q 

as  t  — +  0.  We  can  then  use  the  Maclaurin  series  log(l  -j-  z)  =  z  —  ir^/2  +  z^/3  —  ■  ■  • 
to  obt£Lin 

=  +  (2.3) 

J=l 

as  t  — >  0,  where  the  coefficients,  kj's,  are  the  cumulants  by  definition.  It  is  now 
evident  that  the  cumulant  of  order  r  exists  if  the  moments  of  order  r  and  lower 
exist. 
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Joint  cumulants  that  involve  two  or  more  random  Vtiriables  are  also  of  inter¬ 
est.  Consider  r2indom  variables  such  that  <  oo,  j  ~ 

for  a  positive  integer  n.  Then  all  mixed  absolute  moments  and  cumulants  of 

Xi, . . . ,  Xr,  up  to  order  n  exist.  Let  ^Xi- •  •  • ,  <r)  =  ^jexp  i 

the  joint  characteristic  function  of  X\, . . . ,  Xr,  then  joint  cumulants  of  Xi, . . . ,  Xr, 
cum(  A'l, . . . ,  A'r),  is  given  by  the  coefficient  of  the  •  •  ■  fr/r!  term  in  the  Taylor 
series  expansion  of  log  4>Xi.  Xr(  )  about  the  origin. 

.According  to  the  above  definition,  the  first  four  joint  cumulants  are 
cum(.Vi)  =  £A'i, 

cum(A'i.  A’^i  =  £'[(Xi  -  £.Vi)(A'2  -  £A'2)]  =  Cov(X,.  A'z)^ 
cum(AT.  X2,  X3)  =  E  [(A,  -  £X,)(X2  -  £X2)(X3  -  £X3)] , 

and 

cum(X,.X2.X3,  A'4)  =  £[(A,  -  £X,)(A'2  -  £X2)(X3  -  EX,){X,  -  EX,) 

-  Cov(X,.X2)Cov(X3,X4)  -  Cov(X,,X3)Cov(X2,X4) 
-Cov(Xj.X,)Cov(X2,^3). 

In  general,  the  rth  order  joint  cumulant,  cum(A'’i, . . .  ,Xr),  is  given  by 

cum(Xi . X)  =  ;^(-ir‘(p  -  !)!(£  H  11  (-4) 

jCi'i  ;e*'p 

where  the  summation  extends  over  all  partitions  i/i. . . . ,  i/p,  p  =  1 . r,  of  the  set 

{1. . . .  ,r}. 


Representing  joint  moments  in  terms  of  joint  cumuiants  are  also  useful.  We 


recoid  those  that  will  be  used. 


=  cum(Xi,-Y2)  +  cum(A''i)cum(.Y2), 

E[XiX2X3\  =  cum(Ai,  A''2,X3)  +  cum(Xi)cum(A2,  A3)  +  cum(A2)cum(Ai,  A'’3) 
+  cum(A3)cuin(Ai,  A2)  +  cum(Ai)cum(A2)cum(A3), 
£:[AiA2A3A4]  =  cum(Ai,  A2,A3,  A4)  +  cum(Ai)cum(A2,  A3,  A4) 

4-  cum(A2)cum(Ai,  A3,  A4)  +  cum(A3)cum(Ai,  A2,  A4) 

+  cum(A^4)cum(A'’i,  A'2,  A3)  +  cum(Ai,  A2)cum(A3,  A'’4) 

+  cum( A'’! ,  A3)cum(A^2,  A'’4)  -f  cum(Ai,  A4)cum(A2.  A3) 

4-  cum(A'’i,  A'2)cum(A^3)cum(A4)  4-  cum(Ai,  A3)cum(A2)cum(A4) 
4-  cum(Ai,  A4)cum(A'2)cum(A3)  4- cum(A3,  A4)cum(A'i )cum( A2) 
4-  cum(A2,  A4)cum(A'i)cum{A3)  4-  cum(A2,  A3)cum(Aj)cum(  A4) 
4-  cum(Ai  )cum{A2)cum{A3)cuni(A4). 


For  convenience,  if  i,  j.  and  k  are  positive  integers,  for  random  variables  A', 
y,  and  Z,  we  define 

«:,(A)  =  cum(A, . . . ,  A),  (2.5) 


«OlA,y)  =  cum(A,...,A,y,...,r), 

t  terms  j  terms 


and 


A,  y  z)  =  cum(^. . ,  A,  y, . . . ,  y,  z, . . . ,  z). 


I  terms  j  terms  k  terms 


(2.6) 


We  record  here  some  useful  properties  of  cumulcints;  see  Brillinger  [6],  p.  19. 
Note  that  in  all  cases  p.  r.  and  s  are  positive  integers. 
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(1)  cum(aiXi, . . . ,  arXr)  =  oi  •  •  •  arCum(-Yi, . . . ,  Ar)  for  constants  oi . a,. 

(2)  cum(Ai, . . . ,  A'r)  is  symmetric  in  its  arguments  so  that  cum(Ai, . . . ,  A%)  = 
cum(A^;r, , . . . ,  A„J  for  any  permutation  (tti,  . . . , tt^)  of  (1, . . .  ,r). 

(3)  cum(Ai, . . . ,  Ar)  =  U  if  any  nontrivijil  proper  subset  of  A,’s  are  independent 
of  the  remaining  A,’s.  To  see  this,  we  may  suppose  that  (Ai,...,Ap)  and 
(A'p+i, . . . ,  Ar)  are  independent,  then 

T  P  r 

log  |£  exp  =log{£exp[z^t,A;]}+log{£exp[i  ^  ^,A,]}. 

j=i  j=p+i 

There  will  be  no  {iyti  •  •  •  £  term  in  the  Taylor  series  expansion  on  the  right 
hand  side;  the  desired  result  now  follows  from  equation  (2.4),  the  definition  of 
joint  cumulants. 

(4)  cum(Ai  +  y,  A2, . . . ,  Ar)  =  cum(A'’i,  Aj,. . . ,  A,)  +  cum(y,  A2, . .  • ,  A'r). 

(5)  cum(A'i  +  Oi, . . . ,  Ar  +  flr)  =  cum(Ai, . . . ,  Ar)  when  r  >  2  and  a/s  are  con¬ 
stants.  It  is  sufficient  to  note  from  (4),  that  cum(A'^i  +  ai,  A2, . . . ,  A'r)  = 
cum(Ai,  A2,. . . ,  Ar)  +  cum(ai,  A’2, . . . ,  Ar),  but  the  last  term  is  zero  from 
(3).  The  result  now  follows  from  induction. 
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where  both  second  summations  extend  over  all  nonnegative  integers  i,’s  and 
m’s,  and  positive  integers  jVs,  such  that  it  =  m  and  5Z*=i  f 

(Lukacs  [20],  p.  27). 


(7)  For  a  random  variable  X  the  cumulants  {/Cr(X)}  and  central  moments  {nr{X)) 
satisfy  the  following  relationship. 


Kr  = 


«=0 


‘X- 


'll'" 


(2.10) 


(2.11) 


where  both  second  summations  extend  over  all  nonnegative  integers  *,’s  and 
m’s,  and  integers  jVs,  j,  >  1  for  each  s,  such  that  it  =  ^  and  ]Ct=i  ~ 
r.  To  see  this,  let  X'  =  X  —  EX.  Notice  that  Ki(X')  =  fi\{X')  =  fii{X')  =  0; 
Ki{X')  =  Ki(X)  and  //'(A'')  =  fii{X')  =  fii{X),  for  each  i  >  2.  Applying  (6) 
to  k,(A")’s  and  ^'(vY')’s  yields  the  desired  result. 


(8)  Joint  moments  can  be  representing  by  sums  of  products  of  joint  cumulants. 


(2.12) 

where  the  summation  extends  over  all  partitions  i/i, . . .  ,i/p,  p  =  1, . . . , r,  of 
the  set  {1, . . . ,  r),  and 

K,,,  =  cum(Ao,,...,AE„),  (2,1.3) 

where  the  q/s  are  the  elements  of  u,  (Rosenblatt  [22],  p.  34). 
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2.2.  Formal  Edgeworth  Expansion 

We  now  discuss  some  results  for  asymptotic  expansions  of  distribution  func¬ 
tions.  In  order  to  motivate  the  idea  of  the  expansions,  we  proceed  formally  as 
follows. 

Suppose  Fx,  ^A'l  {«r(A'’)  :  r  >  1},  and  Fy  ,  {Kr(V')  :  r  >  1}  are  dis¬ 

tribution  functions,  characteristic  functions,  and  cumulants  of  random  variables 

X  and  y,  respectively.  Since  $a'(0  =  ^xp  |  Kr( A)(i<)’’/r!|  and  ^y(t}  = 
«r(y)(iOV^!}^  we  have 

^  \r 

*x(t}  =  exp  {  ^[x,(.Y)  -  x.(r))i^}4.K(().  (2.14) 

r=l 

If  the  characteristic  function  is  absolutely  integrable  over  (—00,00),  then  Fy 
is  absolutely  continuous  and 

1  /‘°° 

fy{y)  =  Fy{y)  =  —  J  exp{-ity}^y{t)dt.  (2.15) 

The  density  fy  is  bounded  and  continuous  (Chung  [8],  p.  155).  It  follows  that, 
subject  to  existence, 

D^fyiv)  =  ^  r  e.xp{-j7j/}(-z0^^>>-(0d<.  (2.16) 

where  D  =  d/dy  denotes  the  differential  operator.  Then,  formally,  the  characteris¬ 
tic  function  of  the  distribution  function  ^Fy  will  be  (  —  ity^y{t).  By  this  obser¬ 
vation  and  from  equation  (2.14),  the  uniqueness  property  of  the  Fourier  transform 
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then  yields 

Fx  (i)  =  exp  {  f;ix,(A')  -  (j,),  (2.17) 

r  =  l 

and  similarly 


fxix)  =  exp|^[K,(A')-/Cr(r)]^-^|/y,(x).  (2.18) 

r=l 

The  special  case  where  y  is  a  normal  random  variable  is  most  important.  Let 
y  be  and  set 

a(y)  =  (l/\/2^)exp{-yV2}-  (2.19) 

Then  Y  has  density  function  fyiy)  =  (x~^a((y  —  ^)/<7),  characteristic  function 

$y(f)  =  exp{i/i<  -  and  cumulants  /ci(y)  =  //,  K2{Y)  =  and  Kr{Y)  =  0. 

for  each  r  >  3.  Equation  (2.18)  becomes 

r  /  \  r  ,  «2(A')-(72 

fx  i^)  =  exp  I - - - D  + - - - 

+  (2.2U) 

This  is  the  well  known  Edgeworth  expansion  of  Type  A  (Kendall  and  Stuart  [16], 
p.  170). 

Another  form  of  the  Type  A  series  is  also  of  interest.  First  note  that  the  various 
derivatives  of  fy  in  the  expansion  can  be  expressed  in  the  form  of  Chebyshev- 
Hermite  polynomials.  Specifically,  we  have  (Kendall  and  Stuart  [16],  p.  167) 

(-Dya{x)  =  Hr{x)a{x),  (2.21) 
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where  {Hr{x)  :  i  >  0}  are  the  Chebyshev-Hermite  polynomials  given  by: 

Hoix)  =  1, 

H2{x)  = 

Hi{x)  =  —  3x, 

H^{x)  =  X*  —  6x^  +  3, 

and,  in  general,  the  rth  Chebyshev-Hermite  polynomial  has  the  form 


r^/2i 


H,(x)  =  Y, 


J=0 


{r-2j)\j\2y 


(2.22) 


where  fx]  denotes  the  smallest  integer  which  is  greater  than  or  equal  to  x.  Since 
/y  (x)  =  a((x  —  m)/<7),  from  equation  (2.21), 


{-DYfy  (x)  =  <T-^Hr{{x  -  /i)/^)/^  (x).  (2.23) 

From  equations  (2.18)  and  (2.23)  we  can  see  that  fx  can  now  be  formally  expanded 
in  the  products  of  fy  and  the  Chebyshev-Hermite  polynomials. 

Up  to  this  point  an  underlying  assumption  is  that  the  functions  and 
possess  convergent  Type  A  series.  This  is  not  always  the  case.  For  a  discussion  of 
the  convergence  properties  see  Kendall  and  Stuart  [16],  pp.  173-174,  and  Cramer 
[10],  p.  223. 

For  practical  applications,  however,  it  is  usually  of  little  value  to  know  the  con¬ 
vergence  properties  of  the  expansions.  What  we  are  really  interested  in  is  whether 
a  small  number  of  terms  would  suffice  to  produce  good  approximations  of  the  func¬ 
tions  fx  and  Fx  ■  If  this  is  the  case,  we  would  not  be  too  concerned  about  the 
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convergence  properties.  On  the  other  hand,  if  the  series  actually  converges  but  a 
satisfactory  approximation  can  only  be  obtained  after  a  large  number  of  terms  have 
been  calculated,  then  this  Type  A  series  would  be  of  very  little  use. 

One  important  application  of  the  Type  A  series  is  when  the  cumulants  of  the 
random  variables  have  a  special  structure.  For  example,  suppose  we  are  given  a 
sequence  of  random  variables  {Sn  :  n  >  1},  with  ESn  =  0,  n  =  1, 2, . . . ,  and  for 
each  i  >  2,  the  eth  cumulant  of  5„, 


k-,(5„)  =  0(n),  (2.24) 

as  n  — >  oc.  We  are  interested  in  the  expansion  of  the  distribution  function  of  the 
normalized  random  variable 

Z„  =  (2.25) 

as  n  — >  oo. 

1  /2 

We  use  this  notation:  =  «i(5'n)  and  li  =  «i(Z'n)  =  >  ^2  =  «2(^n)  ~ 

1=0,  and  li  =  «:,(Zti)  =  for  each  i  >  3.  From  equations  (2.24)  and  (2.25), 

it  follows  that,  for  each  i  >  3, 


/.  =  C>(n'-/2).  (2.26) 

Since  Ki(Z„)  =  0  and  K2(Z„)  =  1  we  will  choose  /i  =  0  auid  =  1  in  equation 
(2.18),  the  Edgeworth  expansion  of  Type  A,  and  obtain 

IzM)  =  ^xp  {  -  |d=  +  +  ■  •  •  }a(x).  (2.27) 

Using 

expln'^-^^di  +  n~^d2  +  +  •  •  •} 
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—  1  +  n  di  n  ^  —d\  +  ^2  +  ^1^2  +  <^3  +  ■  •  • ,  (2.28) 

L  J  Lz  J  Lb  J 


and  equation  (2.26),  we  can,  at  least  formally,  approximate  by  the  derivatives  of 

a(-)  up  to  any  order  of  For  our  purpose,  we  neglect  terms  of  order  ©(n”^/^), 

then  we  have,  as  n  — ♦  00 

/z.(i)  =  exp  {  -  +  I^D*  - 

=  Q(i)  +  |//3(z)a(i)  +  a(x) 

0  L  24  ( 2  J 


If  we  apply  similar  argument  to  we  then  have  (Cramer  [10]),  as  n  — >  00 
Fz.(x)  =  4>(i)  -  |/f3(i)a(x)  -  [^*(x)  +  ;|tf5(x)]Q(r) 

Note  that  if  we  only  include  the  first  term  on  the  right  hajid  side  this  is  the  usual 
centred  limit  theorem  while  all  the  remaining  terms,  with  coefficients  of  smaller 
order  of  n,  represent  the  error  terms  in  the  approximation  of  the  central  limit 
theorem. 

In  practical  applications,  we  may  not  know  the  exact  value  for  K2{Sn),  which 
is  the  the  variance  of  5„.  Suppose  that  K2(5’n)  is  an  estimator  for  K2{Sn]  then  a 
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formal  Edgeworth  expansion  of  Z„  =  5„/a;2(5'„) 


1/2 


will  be  of  interest. 


If  /,  =  K:(Zn),  for  each  i  >  1,  satisfies  some  asymptotic  properties,  then  a 
formal  Edgeworth  expajision  can  also  be  obtained.  For  example,  suppose  /j  = 

and  for  each  i  >  3,  /,  =  Then  it  can  be  shown 

that,  as  n  — »  oo  (Kendall  and  Stuart  [16],  p.  176) 


+  ^^3(x)  + 


(.;T)4-„{r-').  (2.31) 


2.3.  Cornish-Fisher  Expansions 

Suppose  that  we  have  a  sequence  of  random  variables  {y„  ;  n  >  1}  which  are 
asymptotically  normal  in  the  sense  that  there  exist  sequences  {/i„  :  n  >  1}  and 
{(7„  :  n  >  1}  such  that  as  n  oo. 


^  N {0,1).  (2.32) 

(^n 


For  large  n,  it  is  often  possible  to  approximate  the  distribution  of  {Yn  —  ^„)/(Tn 
by  a  normal,  but  for  small  to  moderate  n  this  may  not  be  a  good  approximation. 
Under  some  circumstance  we  are  able  to  use  a  polynomial- variate  transformation 
such  as 


ifn  —  ^n.O  "i"  ^n.l  ( 


(Tn 


(2.3u) 
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where  6„,,’s  are  of  order  or  smaller.  By  choosing  6„,,’s  appropriately,  we  may 
make  the  distribution  of  much  closer  to  normal  than  that  of  [Y^  —  Hn)l<^n- 

It  turns  out  that  this  problem  is  embedded  in  the  following  larger  and  more 
natural  question.  Suppose  ^  is  a  A^(0, 1)  random  variable,  Z„  has  distribution  func¬ 
tion  F2^  and  and  its  cumulants  satisfy  equation  (2.26).  If  x  and  ^  are  corresponding 
quantiles  of  and  $  respectively,  namely, 

FzS^)  =  ^(0,  (2.34) 

solving,  at  least  formally,  ^  in  terms  of  x  and  x  in  terms  of  ^  will  be  of  interest. 

We  sketch  how  to  solve  equation  (2.34)  when  can  only  be  estimated. 
The  basic  idea  is  to  approximate  F^^  by  its  Edgeworth  expansion.  Let  z  be  the 
difference  between  distribution  functions  F^^  and  then  as  n  — ►  oo 

=  -[!|H^ix)]a{x)  -  l^H,ix)  +  lH,(i)la(x)  -  ■  ■  (2.35) 

If  we  neglect  terms  of  0(71”*^),  then 

X  =  -l|H,(i)la(x)  -  l^//3(x)  +  ^H^{x)]a{x]  +  o(n-').  (2.36) 

Then  from  equations  (2.34),  (2.36),  and  a  Taylor  expansion, 

(  =  ^-^iF2jx)) 

=  $-'($(x)  +  x) 

=  $-M«»(x))  +  x[<i>-‘r($(x))  +  [$-']"(^»(x)) 
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(2.37) 


To  evaluate  [$-»]'($(•)),  [$-^r($(-)),  and  [^►-'n$(-)),  we 

need  the  following 

observation.  Suppose  ^p{-)  is  a  monotone  function  and  V’(')  =  V’ 

“^(•)  is  its  inverse 

function.  If  the  third  derivatives  of  both  (^(•)  and  ^(•)  exist,  then  it  can  be  shown 

that 

(2.38) 

(2.39) 

(2.40) 

From  equations  (2.21),  (2.38),  (2.39),  and  (2.40),  it  follows  that 

[$"‘]'($(x))  =  l/a(x), 

(2.41) 

[$-‘n$(x))  =  i/^(x)/[a(x)]^ 

(2.42) 

[^-rWx))  =  [3HUx)  -  H2{x)]/[a{x)r. 

(2.43) 

Finally,  observe  that  2  is  a  polynomial  of  again,  we  can  obtain  a  formal 

approximation  of  (2.34)  up  to  any  order  of  Substitute  the  approximation 

(2.36)  into  equation  (2.37)  and  truncate  those  terms  of  order  o(n~^).  We  obtain 

i  =  x-  ^(x^-1)  +  ^(4i^  -  7x)  -  ^(i^  -  3x)  +o(n-^).  (2.44) 

Notice  that  the  three  terms  on  the  right  hand  side  are  of  orders  0(1), 
and  0(n“*),  respect ivel}’. 


18 


In  practice,  it  is  sometimes  more  convenient  to  express  x  in  terms  of  This 
may  be  done  by  a  technique  introduced  by  Cornish  and  Fisher  [9].  First,  let  rj 
represents  the  difference  between  quantile  points  x  and 

r){x)  =  x-i,  (2.45) 

From  equation  (2.44),  it  follows  that 

lii)  =  -  1)  +  -  30  -  1(40  -  70  +  »("-')  =  0(n-'«),(2.46) 

ri'(0  =  |(0  +  j(("  -  1)  -  -  7)  +  o(n--l  =  (2.47) 

I'W  =  I  +  j(0  -  I  +  <-("■■)  =  0(n-''=).  (2.48) 

Then  observe  that  by  a  Taylor  expansion, 
x-,f  =  +  x-<f) 

=  7(0  +  -  07'(0  +  -  0V'(0  -!-••• 

=  7(0  +  {7(0  +  (^-  07(0  +  -  0V(0  +  •  •  •  }»7'U) 

+  jl-iCO  +  (x  -  (WiO  +  5(x  -  0"';"(0  +  •  •  •  }\"(0 
+  (2.49) 

where  the  last  equality  is  obtained  by  substituting  the  value  of  x  —  ^  on  the  right 
hand  side  of  the  second  equality.  Using  this  technique  again  and  again  yields,  at 
least  formally, 

=  7(0  +  7(07'(0  +  +  ^7^(07"(0} 

+  +  ^7'<07'(07"(0  +  l7^(07"'(0}  +■■■■  (2.50) 
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Finally,  substitute  approximations  (2.46),  (2.47)  and  (2.48)  into  (2.50)  and 
neglect  terms  of  o(n~^),  we  then  have 

I  =  {  +  |(£"  -  1)  +  {^(e  -  30  -  ^(2£"  -  5{)}  +  (2.51) 

This  is  the  required  representation  of  i  in  Notice  that  the  three  terms  on  the 
right  hand  side  are  of  orders  0(1),  0(n“*/*),  and  0(n“^),  respectively. 

Both  equations  (2.44)  and  (2.51)  are  approximated  to  order  0(n“^).  Cornish 
and  Fisher  [9]  also  give  higher-order  expansions  and  tables  to  facilitate  their  use. 

The  above  Cornish-Fisher  expansions  (2.44)  and  (2.51)  correspond  to  the 
Edgeworth  expansion  (2.30),  namely,  as  n  —*  oo 


is  useful  when  we  do  not  have  the  exact  value  of  the  variance  of  a  random  variable. 
The  corresponding  Cornish-Fisher  expansions  can  be  obtained  by  applying  the  same 
procedures.  We  list  the  results  as  follows. 


^  =  X- 
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+ 


^(4i^  -  7i)  -  -  3x)  I  +  o(n 


(2.54) 


X  =  ^  + 
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-  3^  -  ^(2^"  -  5^)]  +  o(n-i).  (2.55) 


An  application  of  the  formal  Cornish-Fisher  expansions  (2.44),  (2.51),  (2.54), 
and  (2.55)  is  as  follows.  It  is  quite  often  in  the  simulation  study  that  we  only  know 
a  distribution  function  is  asymptotically  standard  normal  and  we  axe  interested 
in  its  6  quantile  point  xs-  Traditionally  we  are  used  to  let  the  corresponding  quantile 
point  ^4  of  $  as  a  substitute  for  x^,  namely, 

Xs  =  ^6-  (2.56) 

But  according  to  the  Cornish-Fisher  expansions  (2.54)  and  (2.55),  if  we  either 
know  or  have  a  method  to  estimate  the  normalized  cumulants  /,'s  of  Zn,  then  the 
corresponding  (estimated)  adjusted  quantile  points 

/i.(^^)  =  6  +  /7  +  |(^|-l)  (2.57) 

and 


^2(6)  =  6  +  (7  +  ^(^1  -  1)  +  1^6  +  ^((s  -  3^4)  -  ^(2^,'  -  56) 


(2.5S) 


will  be  both  asymptotically  better,  as  estimators  for  if,  than  6- 
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2.4.  Uniqueness  Properties  of  Cornish-Fisher  Expansions 


There  are  two  interesting  interpretations  of  the  Cornish-Fisher  expansions. 
VVe  state  them  below  but  leave  their  proofs  to  the  Appendix. 

Proposition  2.1.  Suppose  ^  and  Xn  are  random  variables,  where  is  N(0, 1) 
and  Xn  satisfies  /ci(A'„)  =  ln,i  =  /C2(A’„)  =  1  -H  /„,2  =  1  +  0(n“^), 

/f3(-Vn)  =  In, 3  =  K4(A„)  =  /„,4  =  0(n“^),  and  /c,(A„)  =  0(n“^),  i  >  5. 

Then. 

( 1 )  there  is  a  polynomial  of  degree  two,  gi(Xn)  =  Xn  +  a„,o  +  such 

that  On.,  =  0(n~^^‘^},  i  =  0,  1,  2,  and  jA:,(0  —  K,(^i(A„))j  =  0(n“^)  for  each  i, 
1  <  i  <  3.  Moreover,  the  coefEcients  an,i's  are  unique  up  to  Specifically, 

fln.o  =  -/n,i  +  (l/6)/„j  +  o(n"^/^),  an, I  =  J 

(2)  If  in  addition,  k,(.Y„)  =  o(n~^),  i  >  5,  then  there  is  a  polynomial  of  degree 
three,  p2(-V„)  =  Xn  +  a„.o  -t-  a„,iA''„  •+  a„,2.\'^  +  a„,3X^,  such  that  an,,  = 

for  each  i,  1  <  i  <  4,  and  |«,(0  ~  ^«(^2(A'„))j  =  c»(n~*)  for  each  i,  1  <  i  <  4. 
Moreover,  the  coeScients  an,i ’s  are  unique  up  to  0(n~^).  Specifically,  a„_o  =  —U.i  + 
(1/6)/„.3  +  o(n-i),  =  -(l/2)/„,2  4  (l/3)/„yn.3  +  (l/8)/„.4  -  {7/36)/^^3  +  o(n-M. 

an,2  =  -(1/6)/„.3  +  o(n~'),  and  a„,3  =  -(l/24)/„,4  +  (1/9)/^, 3  4  o(n“^). 

Proposition  2.2.  Suppose  (  and  are  random  variables,  where  ^  is  :V(0. 1) 
and  Xn  satisfies  ki(A'„)  =  /„,i  =  0{n~^f^),  K2(Xn)  =  1  4  ln,2  =  1  4 
«3(A'n)  =  in, 3  =  0(n"‘/^).  Then, 

(1)  there  is  a  polynomial  of  degree  two,  fii(0  =  ^  +  Qn.o  +  i^n,i(  4  an.2s^-  such 
that  On.,  =  i  —  0,  1.  2,  and  |x,(A„)  —  K,[hi{Q)\  =  0(n"’)  for  each  i. 
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1  <  2  <  3.  Moreover,  the  coefBdents  a„  ,  ’s  are  unique  up  to  SpecificalJy, 

fln.o  =  /n,i  -  (l/6)/n,3  +  a„.i  =  o{n-'^l'^),  and  a„,2  =  (l/6)/„,3  + 

(2)  H  in  addition  K4{Xn)  =  Ua  =  0{n~^),  then  there  is  a  polynomial  of  degree 
three,  /j2((f)  =  +  an,o  +  such  taht  a„,,-  =  0(0“*/^)  for  each 

i,  1  <  i  <  4.,  and  |/c,(A'’„)  —  a:,(/i2(0)I  =  ^or  each  i,  I  <  i  <  4.  Moreover, 

the  coefEcients  a^^i’s  are  unique  up  to  0{n~^).  SpecificalJy,  a„_o  =  L.t  —  (l/C>)/„,3  + 

o{n-'^),  a„,i  =  (1/2)/„,2  -  (l/8)/„,4  +  (5/36)/^,3  +  ofn"'),  a„.2  =  (l/6)/„.3  +  o{n-^), 
and  a„.3  =  (l/24)/n.4  -  (1/18)/^  3  4-  £>(«■*). 

In  Proposition  2.1  we  can  actually  obtain  a  slightly  better  result.  For  example, 
in  part  (1)  we  can  calculate  a„,,’s  so  that  they  will  be  unique  to  order  0{n~^). 
However,  according  to  our  assumption  of  the  cumulants  of  An,  the  neglected  terms 
in  the  formal  Cornish-Fisher  expansion  of  Xn  are  of  order  0(n"').  Thus  calculating 
a,’s  to  a  lower  order  will  not  produce  a  better  approximation,  at  least  orderwise. 
For  this  reason  we  do  not  present  that  result.  A  similar  result  and  observation  also 
apply  to  part  (2). 

One  significant  point  of  Propositions  2.1  and  2.2  *hat  both  the  transforma¬ 
tions  from  a  standard  normal  random  variable  to  a  statistic  and  the  transforma¬ 
tion  from  a  statistic  to  a  standard  normal  random  variable  2ire  unique.  Moreover, 
the  second  and  third  order  polynomial  transformations  coincide  with  the  first  two 
Cornish-Fisher  expansions;  see  equations  (2.54)  and  (2.55). 

.4s  a  consequence  of  an  Edgeworth  expansion,  Cornish-Fisher  expansions  can 
transform  a  statistic  so  that  its  distribution  is  closer  to  a  standard  normal  random 
variable.  Thus  it  can  be  used  to  generate  an  asymptotically  better  confidence 
interval.  Unfortunately,  for  many  statistics  associated  with  real  world  processes,  it 
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is  extremely  difficult  to  prove  the  existence  of  a  rigorous  Edgeworth  expansion.  The 
two  propositions  above,  in  a  sense,  provide  us  an  alternative  approach  for  obtaining 
“formal”  Cornish- Fisher  expansions  for  these  situations. 

2.5.  Cumulants  of  Stationary  Processes 

In  this  section,  we  will  study  some  properties  of  the  cumulants  of  stationciry 
processes. 

Suppose  {A'n  :  n  >  1 }  is  a  discrete  time  strictly  stationary  process  with  mixing 
constants  {a„  :  n  >  1}  such  that  for  any  positive  integers  k  and  n, 

\P{AnB)-P{A)P.iB)\<C‘n 

for  aU  sets  A  and  B,  where  A  €  the  cr-field  generated  by  random 

variables  Ai, . . . ,  Xk,  and  B  €  (T{Xk+ni  -^Jt+n-n?  •  •  •)>  tfie  tr-field  generated  by  ran¬ 
dom  variables  _  If  q„  — »  0  as  n  — ♦  oo,  then  Xk  and  Xk+n  are 

approximately  independent  for  large  n.  In  this  case  the  sequence  {A„  :  n  >  1 }  is 
said  to  be  strongly  mixing  (Rosenblatt  [22],  pp.  63-64). 

Two  discrete  time  strictly  stationary  processes  {A„  ;  w  >  1}  and  :  n  >  1} 
are  jointly  strongly  mixing  with  mixing  constants  {q„  :  n  >  1}  if  for  arbitrary 
constants  Ci  aad  C2,  let  Z„  =  ciA„  +  C2Y„,  for  each  n  >  1,  then  the  process 
{Z„  ;  n  >  1}  is  strongly  mixing  with  mixing  constant  {q„  :  n  >  1}. 

Let  us  define  S„  =  A,  to  be  the  partial  sum,  and  A„  =  S„/n  to  be  the 
sample  mean  associated  with  the  sequence.  The  subscript  n  will  be  omitted  when 
there  is  no  confusion.  We  now  cite  the  following  theorem. 
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Theorem  2.3  (Billingsley  [5],  p.  316).  Suppose  that  the  sequence  {X^  :  n  > 
1}  is  stationary  and  strongly  mixing  with  =  0(n~®)  and  that  EX„  =  0. 

(1)  HEX*  <  oo,  then 

OO 

n  Var[X„]  a"*  =  EXl  +  2  ^  EX^Xk^^ 

k=l 

where  the  series  converges  absolutely. 

(2)  HEX^^  <  oo  and  cr^  >  0,  then  n*l'^X,,l<T  7V(0, 1). 

Corollary  2.4.  Suppose  that  the  sequences  {Xn  :  n  >  1}  and  {Vn  :  n  >  1} 
are  stationary  and  jointly  strongly  mixing  with  q„  =  0(n"®)  and  that  EX^  =  0, 
EYn  =  0,  EX*  <  oo,  and  EY*  <  oo.  Then 

OO  OO 

n  Cov[X„,F„]  ^  cr;fy  =  E[X,Y,]  + 

*=1  *=i 

where  the  series  converges  absolutely. 

.As  noted  in  Billingsley  [5],  p.  316,  in  the  proof  of  Theorem  2.3,  the  conditions 
a„  =  0(n“®)  and  EX^^  <  oo  eire  actually  stronger  than  necessary;  they  are  imposed 
to  avoid  technical  complications  in  the  proof.  One  can  also  see  that  the  sufficient 
condition  for  Vcu:[.Y„]  =  0(n“')  is  weaker  that  that  of  the  central  limit  theorem. 
Actually,  as  we  will  see  in  the  next  theorem,  some  general  properties  about  the 
order  of  the  cumulants  of  the  sample  mean  can  be  obtained. 

Theorem  2.5  (Titus  [25],  p.  16).  Consider  a  stationary  sequence  {.Y„  :  n  >  1} 
such  that:  (1)  for  some  positive  integers  j  and  k,  fjXi  <  oo,  and  (2)  the 
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sequence  is  mixing  mth  a„  =  for  some  e  >  0.  Then  for  constants 

ao,...,ak,  0(n)  as  n  — »  oo.  Consequently,  for  constants 

ao, . . .  ,afe,  Kj{YlLi  J's  0{n^~^)  as  n  oo. 

Instead  of  the  assumptions  in  Theorem  2.5  above,  Titus  [25]  actually  assumes 

that: 


(1 )  For  some  positive  integers  j  and  k,  each  mixed  moment  of  the  form 
E\Xn^  •  ■  •  Xn,  1  is  bounded  for  all  I,  I  <  I  <  4(ji  — 1),  and  (2)  the  sequence 
is  mixing  with  a„  =  for  some  e  >  0. 

For  a  stationary  process,  his  assumptions  and  ours  are  equivcdent.  It  is  easy  to  see 
that  his  assumption  implies  ours.  The  equivalency  can  be  obtained  by  noting  that, 
from  Holder’s  inequality,  •  •  •  Xn,]  <  {E\XnA‘yf‘ ■  ■  •  {E\Xn,\‘y^‘  =  E\XA‘, 

where  the  last  equality  follows  from  stationarity. 

In  both  Theorems  2.3  and  2.5,  the  assumption  that  the  sequence  of  mixing  con¬ 
stants  is  decrecising  faist  enough  is  crucial.  That  assumption  is  true  if,  in  addition, 
the  strictly  stationary  process  {X„  :  n  >  1}  is  either  independent,  moving  average 
of  finite  order,  or  m-dependent.  Moreover,  certain  discrete  time  Markov  chain  has 
mixing  consteints  that  are  decreasing  exponentially  with  time;  see  Billingsley  [5], 
p.  315,  for  details. 

There  are  two  more  results  that  we  will  use  latter  on.  The  first  one  is  a  lemma 
from  Titus  [25]. 

Lemma  2.6  (Titus  [25],  p.  23).  Assume  that  {X'„  :  n  >  1}  is  a  station¬ 
ary  sequence.  Let  p  =  where  each  p,  is  a  positive  integer.  Suppose 
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<  oo,  and  the  sequence  is  mixing  with  a„  =  0{n  ^(p+O)  for  some 
e  >  0.  Then  cum{Sfi\ . . . ,  Sf[>)  is  (9(n“i“(b/2J.p-i+i))  as  n  ^  oo.  Consequently, 
cum(Xl\.  ..,X!)  is  as  n  ^  oo. 

The  second  one  is  a  corollary  of  Theorem  2.5.  Suppose  that  «,(5„)  =  0{n)  as 
n  -♦  oo  for  each  i  <  j,  then  from  properties  (6)  and  (7)  in  Section  2.1,  = 

Otn')  and  /ii(5n)  =  0(nL'/2J)  as  n  oo  for  each  i  <  j.  From  this  observation,  we 
have  the  following. 

Corollary  2.7.  Let  {X„  ;  n  >  1}  he  a  stationary  sequence  such  that:  (1) 
For  some  positive  integers  j  and  k,  <  oo,  (2)  The  sequence  is  rrux- 

ing  with  On  =  for  some  e  >  0.  Then  for  constants  ao,...,ak, 

=  OK)  and  as  n  -  oc.  Con¬ 
sequently,  for  constants  aQ,...,ak,  =  0(1)  and  fl|An)  = 

as  n  — >  oc. 

The  cumulants  of  X„  can  be  calculated  by  the  next  lemma,  which  is  also  from 
Titus  [25]. 

Lemma  2.8  (Titus  [25],  p.  29).  Let  {Xn  ;  n  >  1}  be  a  stationary  process  with 
mixing  constants  {a„  :  n  >  1}. 

(1)  I{E\Xi\‘*  <  oc,  and  a„  =  0(n“‘‘"')  as  n  oc  for  some  e  >  0,  then  as  n  oc, 
«2(X„)  =  n'^C2  +  n-^di  +  o{ti-^),  where 

OO 

C2  =  ^  cum(Xo,X,), 
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oo 

dj  =  -  cum(Xo,X,). 

1=1 

(2)  I{E\Xi\^  <  oo,  and  a„  =  0(n-®-‘)  as  n  —*  oo  for  some  e  >  0,  then  as  n  —*  oo, 

K3(Xn)  =  n"*C3  +  +  o(n“^),  where 

OO 

C3  curn(uV^o> 

i,j=— 00 

00 

^  ^  3i  ^cum(^o»  -^05  X{)  +  cuni(Ao,  Xi,  A^t)| 

t=i  ^ 

00 

-  ^  6(i  +  j)cum(Xo,X„A:j). 

i,j=l 

(3)  If  E\Xi  <  CX3,  and  a„  =  (9(n”*~‘)  asn  00  for  some  e  >  0,  then  as  n  —*  oc, 
K4{Xn)  =  n~^C4  +  n~‘*d4  +  o{n~*),  where 

OO 

C4  =  ^  cum(A'o,X,,X^,.Yjt), 

{,j,l:=-oo 

00 

<^4  =  ~  ^  ^  4z  ^cum(.yo,  .yp,  Aq, .Y,) -|~  cum(Xo,  .y,-,  X,,  .Y,) 

»=i 

00 

-  cum(Xo,Xo,X„Xi) 

1=1 

OO 

-  Y,  12(J  +  j)  [cum(J^o, 

«,j=i 

+  cuiti(-Yq,  Aj,  Aj,  A,.|.j)  -|-  cuni(AQ,  Aq^  Aj,  ) 

00 

-  ^  2^i  +  j  +  k)cum{Xo,Xk,Xj+k,Xi+j+k). 

i,jA=l 
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Corollary  2.9.  Assume  {A„  :  n  >  1}  is  a  stationary  orocess  with  zero  mean  and 
Sn ’s  are  the  pastiaJ  sums: 

(1)  If  K2{Xn)  =  n~^  ■C2  +  n~'^d2  +o(n“^),  then  E{n}l‘^Xn)'^  =  C2  +  n“^ci2  + 

(2)  If  /c,(X„)  =  n-'+i  •  c.  +  n-d.  +  oCn'O  for  2  <  i  <  3,  then  E{n^'^Xr,f  = 

(3)  If  «.(X„)  =  n-+i  •  c.  +  n-'d,  +  o(n-)  for  2  <  i  <  4,  then  E{n^^^Xn)*  = 
3c^  -I-  n~^[6c2d2  +  C4]  +  o(n~^). 

(4)  If  Ki{Xn)  =  •  c,  +  n~‘di  +  o(n~‘)  for  2  <  i  <  5,  then  E(n^^^X„)^  = 

n"^/^[10c2C3]  +  0(n~^^^). 

(5J  If  /Ci(An)  =  •  c,  +  n“'d,-  +  o(n~‘)  for  2  <  i  <  6,  then  E(n^^‘^Xn)^  = 

15c2  +  0(n~^). 


2.6.  Consistency  of  Sample  Moments 

Let  ;  i  >  1}  be  a  discrete  time  strictly  stationary  process.  Assume 
that  {p'  '■  r  >  1}  and  >  1}  are  the  moments  and  central  moments  of 

Xi.  The  {pj.  :  r  >  1}  and  {p,  :  r  >  1}  represent  important  parameters  about 
the  distribution  of  A^i.  Natural  estimators  of  these  parameters  are  given  by  the 
corresponding  sample  moments  of  {X,-  :  i  >  1}.  Thus  may  be  estimated  by 

n 

i=l 

and  pr  may  be  estimated  by 

n 

=  (l/n)^(X,-Xr,  (2.60) 

«=i 
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for  each  r  >  1.  For  convenience,  we  also  denote  =  EXi  and  Xn  =  j. 

In  this  section,  we  follow  the  development  of  Serfling  [24],  pp.  67-71,  to  show  the 
mean  square  consistency  of  these  estimators.  We  begin  w’ith  the  following. 

Proposition  2.10.  Assume  the  sequence  {.V,-  :  i  >  1}  is  stationary  and  strongly 
mixing  with  a,-  =  0{i~^)  and,  for  an  integer  r,  EX*''  <  oo.  Then 
(1)  =//;,• 

(2)  =  0(n-^); 

(3)  MSE{m;,}  =  =  0{n-*); 

(4)  m'„ — ♦  p',.  in  and  in  distribution,  as  n  — +  oo. 

Proof:  (1)  Since  EX'  =  p'  for  each  i,  the  result  follows. 

(2)  Notice  that  the  process  {X'  :  t  >  1}  is  stationary  and  strongly  mixing  with 
o,-  =  0(i~^);  this  result  then  follows  from  Theorem  2.3. 

(3)  This  result  follows  from  (1)  and  (2). 

(4)  This  result  follows  from  (3).  Q 

Preliminary  to  stating  properties  of  the  estimates  it  is  advantageous 

to  consider  the  closely  related  random  variables 

=  (2.61) 

for  each  r  >  1.  Properties  of  m„,r’s  will  be  deduced  from  those  of  the  m*  ,.’s.  The 
same  arguments  employed  in  Proposition  2.10  immediately  yield  the  following. 

Proposition  2.11.  Assume  the  sequence  {X,  :  i  >  1}  is  stationary  and  strongly 
mixing  with  a,  =  C>(i“®)  and,  for  an  integer  r,  EX*'  <  oo.  Then 


(V  ^»^n.r  = 

(2j  Vax{m;,}  =  0(n-M; 

(3)  MSE{m;,,}  =  E{ml,  -  f^r}^  =  0(n-‘); 

(4)  m“ —*  Hr  in  and  in  distribution,  as  n  — ^  oo. 

We  record  the  following  equation  from  Serfling  [24],  p.  69. 

n 

^n.r  =  (l/n)  Y^iXi-JCy 

i=l 

=  (!/")  EE 

•=1  j=0  ^ 

=  E  (2.62) 

j=0 

where  we  define  m*  g  =  1.  Although  analogous  in  form  to  m' ,’s,  the  random 
variables  m„,,.’s  are  much  harder  to  analyze.  Therefore,  instead  of  dealing  with 
mn,r's  directly,  we  exploit  the  relationship  between  m„,r’s  amd  m*  ,.’s.  We  have  the 
following  proposition. 

Proposition  2.12.  Assume  the  sequence  {Xi :  i  >  1}  is  stationary  and  strongly 
mixing  with  =  0(t~®)  and,  for  an  integer  r,  EX*'"  <  oc.  If,  in  addition, 
E{m\y  =  for  each  j,  0  <  j  <r,  then 

(1)  Em„,r  =  Ht  +  0{n~^); 

(2)  Var{m„,r}  =  0(n-*); 

(3)  MSE{m„,r}  =  =  0(n"*); 

(4)  mn,r  Hr  i^  aod  in  distribution,  as  n  oo. 
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Proof:  (1)  Utilize  equation  (2.62)  to  write 


-  Mr 


SC)'-'- 


Now,  observe  that 

=  ('/"“)£{  -  /')'-■  -  >‘)} 

.=i  i=i 

oo 

=  (l/n)  Y,  E{{Xo-fiy-\Xi-fi)}  +  0{n-^) 

|  =  — OO 

=  0(n-^), 

where  the  second  equality  is  obtained  from  the  fact  that  {(X^  —  m)  •  *  —  0  zero 
mean,  {(Xi  -  :  i  >  1}  and  {(X,-  -  m)  :  *  >  1}  are  jointly  strongly  mixing,  the 

fourth  central  moments  of  both  (Xj  —  and  (Xi  —  m)  are  both  finite,  and  then 
by  an  application  of  Corollary  2.4.  For  each  j  >  2,  use  Holder  s  inequality 


By  application  of  Minkowski’s  inequality 


=  [E{Kl/n)  £«-/')'■']} 

tsl 

<(l/n)£[£:(|Jf,-/ir)] 

=  0(1), 


.  1 '•/(r-j)T  (r-j)/r 


[r-j)/r 


<=l 


and  the  assumption  that  E{Tn*  iy  =  we  obtain 


Tj/’’ 
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=  0(n-»). 

The  results  then  foUow. 

(2)  Writing  Var {m„,r}  =  we  seek  to  compute  E{m‘^^}  and 

combine  the  result  in  (1).  To  do  this,  we  need  to  compute  quantities  of  the  form 

i2,  0  <  ji,  jj  <  r.  For  ji  =  72  =  0,  we  have  =  /C2{ni;.r}  +  But 

«2{m;_r} 

1  =  1 

oo 

=(l/n)  Y.  £{[(A:o-/<)'-/^r][(Ji:^-/i)'-/i']}  +  0(n-") 

i=— 00 

=0{n-% 


where  the  equalities  are  justified  by  the  moment  and  mixing  conditions  as  well  as 
an  application  of  Corollary  2.4.  For  (jx, 72)  =  (0, 1)  or  (1,0),  we  have 

^{^n,T-iK,TK.i}  =  cum(m;_,_j,m;,„m;i) 

+  cum(m;,,._i)cum(m;,^,m;_i) 

+  cum(m;  Jcum(m;^_i,m;_i) 

+  cum(m;,i)cum(m;,_i,m;_J 
+  cum(m;,,_i)cum(m;  Jcum(m;_i). 


The  first  term  is  at  most  0{n~^)  by  Theorem  2.5.  Each  of  the  next  two  terms  is  at 
most  0(n”*),  which  can  be  proved  in  a  similar  fashion  as  in  (1).  The  last  two  terms 
on  the  right  hand  side  are  both  zero  since  Em*^^  =  0.  Thus  £’{m*  ,.m*  j }  = 
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0(n  ^).  Finally,  for  ji  +  >  2,  by  a  similar  application  cf  the  Holder’s  inequality 

and  the  Minkowski's  inequality,  we  have 


The  desired  result  now  follows. 

(3)  This  result  follows  from  (1)  tind  (2). 

(4)  This  result  follows  from  (3). 


=  Oin-^). 


□ 
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Chapter  3 


Johnson- Glynn  Pivots  for  Stationary  Processes 


In  order  to  obtain  a  more  accurate  confidence  interval  it  is  necessary  that  we 
first  obtain  the  formal  Edgeworth  and  Cornish- Fisher  expansions  of  the  saimple 
statistic  of  interest.  However,  to  obtain  these  expansions  we  have  to  calculate  and 
estimate  vajious  moments  and  cumulants  of  the  sample  statistic.  In  this  chapter, 
we  will  use  this  approach  to  obtain  the  first  and  second  order  Johnson-Glynn  pivots 
and  the  associated  confidence  intervals  for  the  batch  means  method. 

The  organization  of  this  chapter  is  as  follows.  In  Section  3.1,  we  will  review  the 
traditional  batch  means  method.  Section  3.2  shows  the  calculation  of  some  mixed 
cumulants  related  to  the  batch  means  method.  The  cumulants  of  the  sample  t- 
statistic  and  the  first  and  second  order  Johnson-Glynn  pivots  for  the  batch  means 
method  will  be  presented  in  Section  3.3.  Section  3.4  derives  two  new  confidence 
intervals  from  the  Johnson-Glynn  pivots  and  shows  that  the  increase  of  the  length 
of  the  confidence  intervals  are  asj'mptotically  negligible.  Section  3.5  compares  the 
computational  efficiency  of  the  traditional  batch  means  method  with  that  of  the 
first  and  second  order  Johnson-Glynn  pivots. 

3.1.  Batch  Means  Method 

Suppose  {X,  :  i  >  1}  is  a  discrete  time  strictly  stationwy  process.  We  are 
interested  in  obtaining  a  point  estimate  and  a  confidence  interval  for  the  stationary 
mean  EX ,  where  X  heis  the  same  distribution  as  Xi.  In  this  section,  we  will  review 
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the  traditional  batch  means  method.  Without  loss  of  generality,  we  assume  that 
:  i  >  1}  has  zero- mean. 

To  show  that  a  confidence  interval  for  EX  can  be  constructed,  we  follow 
the  development  of  BriUinger  [6].  However,  our  assumptions  about  the  stationary 
processes  are  slightly  different.  Specifically,  BriUinger  {[6],  p.  419)  assumes  that: 

The  discrete  time  process  {A,  :  i  >  1}  is  stationary,  continuous  in  mean, 
and  for  each  k,i  >  1,  and  ji,  ■ .  ■  ,jk  >  0,  its  cumulant  of  order  k  +  l  exists 
and  satisfies  (cum{A,,  A, . . ,  A,+j^}|  <  I*,(l  •  •  •  (1  +Jk)~^>  ^or 

some  finite  positive  constant  Lk. 

In  our  discussion,  instead  of  the  assumption  that  the  cumulant  functions  decrease 
reasonably  fast  as  the  arguments  become  far  apart,  we  assume  that  the  station¬ 
ary  processes  satisfy  certain  moment  conditions  and  are  strongly  mixing  with  the 
mixing  constamts  q,’s  decrease  as  i  — >  oc.  Notice  that  both  assumptions  require 
that  sufficient  number  of  the  moments  of  the  random  variable  A  are  finite  and  also 
require  that  the  values  of  the  stationary  process  at  a  distance  from  each  other  are 
only  weakly  statistically  dependent.  Both  assumptions  are  true  if,  in  addition  to 
their  respective  moment  conditions,  the  strictly  stationary  process  {A'’,-  :  i  >  1}  is 
either  independent,  moving  average  of  finite  order,  or  m-dependent. 

Notice  that  in  the  following  discussion,  only  discrete  time  strictly  stationary 
processes  are  considered;  however,  the  results  for  continuous  time  processes  can  be 
obtained  in  a  similar  fashion  (see  BriUinger  [6]). 

Let 

m 

A,„  =  (l/m)^A,  (3.1) 

«=i 
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be  the  sample  mean.  From  Theorem  2.3,  we  know  thai  if  the  mixing  constants  q, 
are  0{i-^)  as  i  ^  oo,  EX}^  <  oo,  and  =  EXf  +  2  EX^Xi+i  >  0,  then 
Xmlicr Ivn}!'^)  converges  weakly  to  -V(0, 1).  Thus  if  we  can  estimate  a  then  we  can 
use  Xm.  as  a  point  estimate  and  construct  a  confidence  interval  for  EX. 

Suppose  the  time  interval  (l,m]  is  split  into  n  intervals  of  length  b  =  [m/nj 
each  and  denotes  the  mean  for  the  ith  interval,  i.e. 

i',  =  i  (3.2) 

j=(.-l)6+l 

for  2  =  1, ....  n.  For  the  Batch  Means  Method,  n  is  the  number  of  batches,  6  is  the 
batch  size,  and  the  Yi's  are  batch  means.  The  basic  id^a  for  batching  is  as  follows. 
We  have  a  stationary  process  whose  values  at  a  distance  axe  only  weakly  statistically 
dependent.  By  batching,  we  transform  the  original  data  sequence  to  another  one 
and,  hopefully,  the  dependency  of  that  sequence  will  decay  faster.  As  a  matter  of 
fact,  we  can  show  that  the  batch  means  are  asymptotically  independent,  in  the  sense 
that  \P{A  r\  B)  —  P{A)P{B)\  — >  0  for  all  sets  A  and  B,  where  A  G  (7(yi,. . .  ,Yk), 
and  B  G  cr{Yk+n, . . .),  as  n  — >  oo.  We  state  the  following  proposition,  whose  proof 
will  be  given  at  the  Appendix. 

Proposition  3.1.  Suppose  {Ai  :  i  >  1}  is  a  discrete  time,  strictly  stationary 
stochastic  process.  For  a  fixed  positive  integer  b  and  a  real  function  g:  TZ^  — >  'R., 
where  TZ  is  the  set  of  real  numbers,  let  Y{  =  A(,_i)6+i,  A(,_j)6+25  •  •  •  i 

i  >  1.  Then 

(1)  The  process  {V]  ;  i  >  1}  is  strictly  stationary. 

(2)  Suppose,  in  addition  to  being  stationary,  {X,  :  i  >  1}  is  strongly  mixing  with 
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mixing  constants  {q,  :  i  >  1}.  Without  loss  of  generality,  we  can  assume  that 
the  sequence  of  mixing  constants  {q^  :  i  >  1}  is  a  nonnegative  and  nonincreasing 
sequence  of  i.  Then  {Vj  :  z  >  1}  is  strongly  mixing  with  mixing  constants  {^i  :  i  > 
1},  where  for  each  i  >1, 

(3)  For  an  increasing  function  h-.'lV'  TV',  where  TV  is  the  set  of  nonnegative 
real  numbers,  h(ai)  <  oo  if  and  only  if^l^i  <  oo- 

(4)  Suppose,  in  particular,  Yi  =  (1/h)  -Y(,_i)6+j  for  each  z  >  1.  Then  E\Xi  p  < 

oo  implies  that  T’|Vip  <  oc. 

Basically,  Proposition  3.1  states  tnat  the  stationarity,  mixing,  and  moment 
conditions  of  a  stationary  process,  as  well  as  the  finite  summability  of  the  mixing 
constants  are  preserved  after  batching.  We  then  have  the  following  result. 

Theorem  3.2  (Brillinger  [6],  p.  420).  Suppose  {Xj  :  z  >  1}  is  stationary  and 
strongly  mixing  with  a;  =  0(i“®)  and  that  EXi  =  0  and  EX}'^  <  oc.  Then  for 
fixed  n,  the  Vi’s  are  asymptotically  independent  and  individually  asymptotically 
normal  with  mean  EX  and  variance 


(l/b)<r^  =  {\lb){EXl4-2Y,EX^X,^,), 

fc=i 


as  m  oo. 


Since  the  Y,’s  are  tisj-mptotically  independent  and  and  individually  asymptot¬ 
ically  normal,  the  next  natural  question  is  whether  we  can  use  the  usual  approach 
in  analyzing  i.i.d.  normal  random  samples  to  generate  a  confidence  interval  for  )  . 
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which  is  Xm-  Define 


=  (3-3) 

^  •=! 

as  the  sample  variance  for  {Yi}  and  let 

.  _Xm- EX 

We  have  the  following. 

Corollary  3.3  (Brillinger  [6],  p.  421).  Suppose  {X,  ;  i  >  1}  is  stationary  and 
strongly  mixing  with  a,  =  0(z“®)  and  that  EXi  =  0  and  EXp  <  oo.  Then  for 
fixed  n,  the  random  variable  t’  j,  converges  to  a  Student-t  random  variable  with 
n  —  1  degrees  of  freedom  in  the  limit  as  m  —*  oo. 

Since  a  Student-t  random  variable  with  n  —  1  degrees  of  freedom  converges  to 
A'^(0, 1)  as  n  ^  oc,  we  can  use  the  quantile  points  of  the  Student-t  random  variable 
or  the  standard  normal  random  variable  to  construct  a  confidence  interval  for  EX. 

3.2.  Cumulants  in  Batch  Means  Method 

The  basic  idea  of  the  batch  means  method  is  to  choose  the  batch  size  large 
enough  so  that,  hopefully,  the  correlation  between  the  batch  means  becomes  small 
and  each  batch  means  is  reasonably  close  to  a  normal.  Then  the  confidence  interval 
methods  for  the  i.i.d.  normal  random  variables  cein  be  appUed  (see  Corollary  3.3). 
Thus  the  properties  such  as  joint  cumulants  and  moments  or  the  mixing  constants 
of  the  batch  means  are  important  in  the  batch  means  method.  In  this  section, 
various  cumulants  related  to  the  batch  means  will  be  given. 


(3.4) 
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Recall  that  for  each  j,  1  <  j  <  n,  X(,  -i)6+j  is  the  nh  batch 


mean.  We  also  define 


<Ti=m-  Var(A'm), 


\  =J^_7 

"""K/i) 


We  will  compute  the  order  of  various  mixed  cumulants  of  X„,,  (1/n) 

and  ^n,b-  Id  the  following  derivation,  we  w'iU  assume  that  the  process  {X,  :  1  <  i  < 
m)  satisfies  sufficient  mixing  and  moment  conditions  so  that  there  wiU  not  be  any 
problem  regarding  convergence  of  infinite  sums.  These  conditions  can  be  obtained 
by  applying  Theorem  2.5.  Below  we  list  the  asymptotic  order  of  the  cumulants, 
while  relegating  the  derivation  to  the  Appendix.  Each  of  the  results  below  requires 
certain  order  conditions  on  mixing  constants  and  moments.  As  these  are  quite 
complicated,  they  will  only  be  given  in  the  Appendix. 

Cumulants  of 


K2{Xr.)  =  0(m-») 
=  0(m"^) 
Ki{Xm)  =  0{m-^) 

K,(Xm)  =  0{m-*) 


Cumulants  of  X, 


=  0[m 
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«2(X^)  =  0(m-2) 

Ksixl)  =  0{m-^) 

K,(xi)  =  0{m-*) 

Mixed  Cumulants  of  Xm  and  X^ 

=  0(m-2) 

»,,!(X„,X^)  =  0(m-=) 

=  0(m-^] 
'‘2.2(X^XJ  =  0(m-“) 
»uX,Xl)  =  0(m-*) 
-c<,,(X„,Xl)  =  0(m-*) 
«3,2(X„X)  =  0(™‘‘) 
»2,3(X„,Xi)  =  0(m-<) 
K„(X„,Xi)  =  0(m-‘) 
««(X„,XL)  =  0(m-‘) 
«33(X„X)  =  0(m-‘) 
«.,,(X„,XiL)  =  0{m-«) 

«5,2(X„X)  =  0(m-«) 

'C7,i(X„,Xl)  =  0(m-') 

Cumulants  of  (1/n) 

«.((!/-.)  E”=.^’)  =  0(6-') 
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«=((!/»)  E”=.v;^)  =  0(n-')0(i-“) 

«-((!/")  i;"«l?)=0(n-=)0(t-‘) 

Mixed  cumulants  of  Xm  and  (1/n)  Yl'i=i  K* 

x.,i(X„,  (1/n)  Y?)  =  0(n-')0(6-=) 

»2,i(X„,  (1/n)  V;“)  =  0(n-“)0(4-^) 

n.nCX™,  (1/n)  £?=.  =  0(n-=)0(i.-=) 

n3..(X„,  (1/n)  n.,  Vn  =  0(n-3)0(i-S) 
X3,3(X„,  (1/n)  E”,.  Y?)  =  0(n-=)0((.->) 
«..3(X„,  (1/n)  V(^^)  =  0(n-»)0(4-) 

-=,,l(X„,  (1/n)  el,  =  0(n-‘)0(6-) 

X3,3(X„,(l/n)  el,  1^’)  =  0(n-‘)0(6-) 
X3,3(X„.{l/n)  el,  1?)  =  0(n-*)0(6-<) 

x„(X„,(l/n)  el,  y?)  =  0(n-')0(6-=) 

»,.3(X„,  (1/n)  el,  5?)  =  0(n->)0(i.->) 
»3.3(X„,  (1/n)  el,  y?)  =  0(n->)0(4-») 

«e.,(X™,  (1/n)  el,  =  O(n-«)0(i,-®) 
«,,3(X„,  (1/n)  el,  1^’)  =  0(n-«)0(6-') 
«T,,(X„,  (1/n)  el,  i?)  =  0(n-’)0(4-T 

Mixed  cumulants  of  and  (1/n)  EL, 

«,,,(Xi,(l/n)  el,  y?)  =  0(n-^)0(6-=) 
n3..(X^.(l/n)  el,  y?)  =  0(n-’)0(4-=) 
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«i,2(Ji:L,(l/n)  vn  =  0(n-=)0(6->)  +  0(n->)0(6-<) 


Mixed  cumulants  of  -Ym,  A'^,  and  (1/n)  2?=! 

au..(X„,Ai,(l/n)Er=ii^^)  =  0(n-»)0(6-») 

X2,..i(7«.XL(l/n)  E”=.  '?)  =  0{->-^)0(4-=) 

x..2,i(Xn>,Xl,  (1/n)  el,  yn  =  0(n-“)0((.-) 

y^)  =  o(n-’)0(i.-‘) 

X3,.,i(7n..XL  (1/n)  Er=i  yn  =  0(n-‘)0(i.--) 
x,,!..(X„,XL,  (l/n)  EL,  5?)  =  0(n-)0(6-*) 

X2,,,2(X„,XL  (1/n)  el,  1^")  =  0(n-')0(6-) 
n.,,..(X„,Xl,(l/n)  EL,i^^)  =  0(n-=)0(6-“) 
«3,!,,(X„,Xi,  (1/n)  el,  y^  =  Cl(n->)0(6-'>) 
«3,,,2(X„.Xl,  (1/n)  el,  If)  =  0(n-»)0(4->) 


Cumulants  of  A„,{, 


«i(A„,i)  =  0(6-1) +  0(n-') 

«2(A„,i)  =  O(n-i) 

/C3(A„,6)  =  0(n-*) 


Mixed  cumulants  of  A„,6  and 


/cu(A„,6,X„,)  =  0(n-i)0(6-i) 
K2,i(A„,6,X^)  =  0(n-2)0(6-i) 
fCi.2(A„.6,X^)  =  0(n-2)0(6-i) 
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/£3.x(A„.,,X^)  =  0(n-3)0(6-^) 
«2.2(A„,i,X„.)  =  0(n-3)0(6-i) 
«i.3(A„,fc,X^)  =  0(n-3)0(6-2) 
K3.2(A„.j,X^)  -  C»(n-‘)0(i-») 
«2.3(A„.6,X„)  =  0(n-‘‘)0(6-=*) 
«i,4(A„,6,X^)  =  (9(n-‘)0(i-3) 
/C3,3(A„,i,X^)  =  0(n-5)0(6-2) 
K2,4(A„,t,X„)  =  0(n-5)0(6-3) 


3.3.  Johnson-Glynn  Pivots 


Suppose  a  discrete  time  strictly  stationary  process  {Xi  ;  i  >  1}  is  strongly 
mixing  with  mixing  constants  {q,-  :  i  >  1}.  We  are  interested  in  applying  the  ideas 
of  the  Johnson-Glynn  pivots  to  obtain  a  confidence  interval  for  the  stationary  mean 
EX,  where  X  has  the  same  distribution  as  X,.  Again,  we  assume  that  {X,  :  i  >  1} 
is  zero-mean. 

Our  starting  point  is  the  batch  means  method,  which  we  intend  to  modify 
to  obtain  new  confidence  interval  methods.  The  basic  idea  about  batch  means 
method  has  been  discussed  in  Section  3.1.  Recall  that,  for  each  i,  1  <  i  <  n, 

y<  =  K.,»  =  (l/n)EL.(5^  <  =  m  .  Vai(X„),  and 

=  Ki,i/(a^/6)  —  1;  see  equations  (3.2),  (3.5),  (3.6),  and  (3.7).  We  also  define 


t 


n,b  — 


Xm-EX 

(V;.6/n)*/2’ 


(3.8) 


In  this  section,  we  will  demonstrate  how  to  calculate  the  first  four  cumulants  of 
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^n,b‘ 

To  handle  various  cumulants  of  tn,h,  note  that,  from  equation  (3.7),  = 

{a‘^/b){l  -h  From  Lemma  A.  12  in  the  Appendix,  A„^  converges  to  0  in  L^, 

as  n,6  — »  oo.  This  implies  that  A„,(,  is  small  for  large  n  and  6.  Simple  algebra  and 
a  Taylor  expansion  then  yields 

y—  =  ~  +  ^n,i>  -  ^n,6  +  ■  •  •  (3.9) 

SO  that,  from  equation  (3.8), 

=  — (X,„  -  £X){l  -  iA„.,  +  -A^,,  -  +  •••}.  (3.10) 

It  is  clear  that  if  we  can  calculate  the  various  mixed  cumulants  of  (X^  —  EX), 
(X„i  — £’,Y)A„,6,  {Xm—EX)Alf,, ,  etc.,  then  we  have  obtained  various  cumulants 
of 

The  whole  procedure  is  very  technical  and  the  exact  results  are  very  involved. 
Basically,  the  idea  is  as  follows.  Notice  that  we  assume  that  :  1  <  i  <  m} 
is  zero-mean.  We  need  to  calculate  various  mixed  cumulants  of  Xm,  XmA„,6, 
Xm  A^  Xm  A^  j,  . . . ,  etc.  But  those  cumulants  can  be  obtained  after  calculating 
various  mixed  cumulants  of  Xm  and  A„,4.  Next,  observe  that  from  equation  (3.7). 
A„,4  =  [Ki,6/(a^/6)]  -  1,  but  from  equation  (3.5),  V;,6  =  (l/«)  = 

(l/n)  —  X^,  so  that  we  have  to  C2dculate  mixed  cumulants  of  X^,  X,„,  and 

(1/n)  which  in  turn  can  be  solved  by  first  calculating  mixed  cumulants 
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of  Xm  (!/”)  Z^"=i  Notice  that  the  order  of  each  of  these  cumulants  and 
mixed  cumulants  mentioned  above  have  been  discussed  in  Section  3.2.  Applying 
these  ideaSj  we  derive  the  following. 

Proposition  3.4.  Assume  that  {Xi :  t  >  1)  is  a  discrete  time,  strictly  stationary 
stochastic  process  with  zero  mean.  Then 

(1 )  If  £'|Xi|*  <  oc,  and  the  sequence  is  mixing  with  a„  =  0(n~^),  then  «i(t„,(,)  = 

-(l/2)(  (1/n)  ZL, 

(2)  If  E\X\  1*^  <  oc,  and  the  sequence  is  mixing  with  an  =  0{n  *),  then  Ksltn^b)  = 
1  +  (blcifK,((lln)  ZU  y?)  -  {mblai)K,,,(Xn.,(\ln)  ZU 

-  mh)  EL,  W/4)  -  11  +  3/n  +  0(r.-')  +  0(6-). 

(3)  If  £|Xi  1^®  <  cc,  end  the  sequence  is  mixing  with  On  =  0(n"**),  tiien  K2{in,b)  = 

WT^MH'ts(Xn,)-3{ci/m)K„(Xn.Mn)ZUy?n‘^im+o(n-''Mb-''^)- 

(4)  If  E\Xi\'^°  <  oo,  and  the  sequence  is  mixing  with  Qn  =  then  K4(tn,6)  = 

-6(m6M)«2a(^m,  (1/n)  E?=i  i?)  +  3(6/0^ )2K2((l/n)  +  12/n  + 

It  can  also  be  shown  that  we  have  «i(t„^j)  =  K2(^n,6)  = 

1  +  0(n~^)  +  0(6"^),  K3(t„.6)  =  0(n~^/^)0(b~^/^),  and  K4(i„,i,)  =  0(n~'^).  We  also 
note  that  the  0(6'*)  term  in  would  disappe^lr  should  {1^  :  1  <  i  <  m}  be 

independent. 

The  first  four  cumulants  of  a  standard  normal  random  variable  are  0,  1,  0, 
and  0.  K  we  want  to  use  a  normal  approximation  to  generate  a  confidence  interval, 
we  would  like  to  have  the  cumulants  of  t„,(,  to  be  as  close  to  those  of  the  standard 
normal  random  variable  as  possible. 

This  idea  follows  from  Edgeworth  expansion  in  the  following  sense.  In  the 
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Edgeworth  expansion  for  the  distribution  function  of  t„,i„  the  error  terms  have 
coefficients  in  terms  of  the  cumulants  of  namely,  K2{tn,b)  —  1,  Ksftn.fc). 

and  see  equation  (2.31).  Thus  if  «2(fn,6)  —  1,  and  K4(t„,6) 

are  close  to  zero,  the  effects  of  the  error  terms  will  be  small.  Notice  that  the 
first  four  cumulcints  of  a  stcindard  normal  random  viiriable  are  0,  1,  0,  and  0.  We 
remark  parenthetically  that  requiring  Ki{tn,b)  —>■  0,  K2{tn,b)  —  1  — ►  0,  >  0. 

and  /C4(f„,i,)  — >  0  is  the  s£ime  as  requiring  that  the  first  four  cumulants  of  are 
close  to  those  of  a  standard  normal  random  variable. 


The  above  argument  leads  to  two  naive,  but  natural,  approaches  for  the 
choice  of  the  batch  size:  (1)  Minimizing  maxi<,<4  |k, •(<„,!,)  —  (2)  Minimizing 

Z!,=i[«>(^n.(>)  -  ««(s)]^-  From  Proposition  3.4,  Ki(tn,b)  -  «i(0  = 

«2(fn.i)  -  «2(<f )  =  0(n"^)  +  K3(tn,b)  “  ^3(0  =  and  K4{t„,b)  -  K4U)  = 

0(n"').  In  botn  cas'“s  we  can  see  that  the  optimal  b  and  n  should  be  chosen  such 
that  6  and  n  are  of  the  same  order,  namely,  both  b  ~  and  n  ~ 

as  m  — >  00.  This  relationship  of  the  batch  size  and  the  number  of  batches  gives  us 
a  tn,b  which  is  closer  to  a  standard  normal  random  variable  in  the  sense  that  the 
differences  between  their  respective  first  four  cumulants  are  smaller.  Hence  we  feel 
this  is  a  preferred  choice  of  batch  size  for  the  traditional  batch  means  method. 


•A-fter  calculating  various  cumulants,  the  formal  Edgeworth  expansion  and 
Cornish-Fisher  expansion  can  be  obtained.  Notice  that  for  fixed  m,  the  choice 
of  6  and  n  will  affect  both  formal  expansions.  For  n  ~  0{m}l'^)  and  b  ~ 
the  formal  Edgeworth  expansion  is 


-  «l(<n,6))<A(x)  +  ( 


^4i^n,b)  ^2i^n,b)  1 

8  2 


)x<p{x) 


x^o{x) 


24 


x^<f>{x)  +  o(n"'). 


(3.11) 
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Again  we  can  have  the  following  new  uniqueness  properties  of  the  Cornish- 
Fisher  expansions  for  the  batch  means  method.  The  proofs  of  the  following  two 
propositions  axe  similar  to  that  of  Propositions  2.1  and  2.2  iind  thus  omitted. 

Proposition  3.5.  Suppose  ^  and  are  random  variables,  where  ^  is  A’(0,1) 
and  X„  satisfies  Ki(Ar„)  =  /„,!  =  K2(A’„)  =  1  +  /n,2  =  1  + 

/i3{Xn)  =  ln,3  =  «4(.Y„)  =  /„.4  =  0(n-^),  and  K,(X„)  =  0(n-^),  i  >  5. 

Then, 

(1)  there  is  a  polynomial  of  degree  two,  gi{Xn)  =  Xn  +  a„,o  +  +  an,2-Y^, 

such  that  an,i  =  i  =  0,  1,  2,  and  |«i(0  “  «»(5i(-^n))|  =  0{n~^)  for 

each  i,  1  <  *  <  3.  Moreover,  the  coefEcients  an.,’5  are  unique  up  to  0{n~^l^). 
Specifically,  a„,o  =  -/„,i  +  (l/6)/„,3  +  o(n~^/^),  a„,i  =  -(l/2)/„,2  +  o{n~^^'^),  and 
an.2  =  -(1/6)/„,3  +  o(n~'/2). 

(2)  If  in  addition,  /c,(A'„)  =  o(n“'),  i  >  5,  then  there  is  a  polynomial  of  degree  three, 

52(A'„)  =  +  a„,o  +  such  that  a„,,-  =  (9(n-^/2)  for  each 

i,  1  <  j  <  4,  and  |K,((f)  -  «,(52(A''„))|  =  o(n“*)  for  each  i,  1  <  i  <  4.  Moreover,  the 
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coefEcients  a„,i’s  are  unique  up  to  0(n~^).  SpeciEcaJly,  a-,i,o  =  -/«,!  +  (l/6)/„,3  + 
o(n~'),  On, I  =  — (1/2)/„,2  +  (l/8)/n,4  +  On, 2  =  “ (l/6)/„^  +  o(n~'),  and 

On,3  =  ■~(1/24)/„,4  +  o(n“'). 

Proposition  3.6.  Suppose  ^  and  X„  are  random  variables,  where  ^  is  :V(0,1) 
and  Xn  satisfies  «i(X„)  =  /„,i  =  0{n~^f^),  K2{Xn}  =  1  +  ln,2  =  1  +  0(n“'''2), 

«3(X„)  =  In, 3  =  0(n-l/2) 

(1)  there  is  a  polynomial  of  degree  two,  hi{()  =  i  +  a^.o  +  Qn.ii  +  such  that 

(in, I  =  i  =  0,  1,  2.  and  !«,(.Y„)  — =  0{n~^)  for  each  i,  1  <  ?  <  3. 

Moreover,  the  coefficients  a„^’s  are  unique  up  to  0(n~^^^).  Specifically,  a„,o  = 
/„,i-(l/6)/„,3  +  o(n~^/^),  a„.i  =  (l/2)/„,2+o(n~*/2),  and  a„,2  =  (l/6)/„,3  +  o(n“^/2). 

(2)  If  in  addition  K4(Xn)  =  ln,4  =  0(n~^),  then  there  is  a  polynomial  of  degree 

three,  /i2(0  =  ^  +  a„,o  +  Qn.i^  +  an.2^^  +  an,3i^,  such  taht  a„,i  =  0(n~^^^)  for 
each  1,  1  <  i  <  4,  and  j/c,(A’n)  —  Ki(/i2(0)l  ==  h  1  ^  <  4. 

Moreover,  the  coefficients  a„^,’s  are  unique  up  to  0(n~^).  Specifically,  a„,o  =  U,i  — 
(1/6)/„.3  +  o(n"^),  a„,i  =  (l/2)/„.2  -  (l/8)/„,4  +  an.2  =  (l/6)/„,3  +  o(n~^). 

and  a„,3  =  (l/24)/„,4  +  o(n-^). 

3.4.  Confidence  Intervals  Generated  from  Johnson-Glynn  Pivots 

In  this  section,  we  will  generate  new  confidence  interval  methods  from  the 
Johnson-Glynn  pivots  and  then  show  that  the  increetses  of  the  lengths  of  the  new 
confidence  intervals  generated  by  the  first  and  second  pivots  are  asymptotically 
negligible. 


First,  we  discuss  how  to  use  the  formal  Cornish- Fisher  expansions  (3.12)  and 
(3.13)  to  produce  confidence  intervals.  In  our  discussion,  we  will  assume  that  both 
b  ~  0{vn}^^)  and  n  ~  0(m^/^),  namely,  b  and  n  approach  to  infinity  roughly  at 
the  same  rate.  We  remaxk  that  our  results  in  previous  two  sections  are  stated  in 
orders  of  n  and  b  and  allow  n  and  b  to  vary  individually.  The  current  choice  for 
b  ~  0{m}l^)  cind  n  ~  0(m^/^)  here  illustrates  how  confidence  intervals  can  be 
generated.  Other  choices  of  n  and  b  can  be  be  treated  in  a  similar  fashion. 

Notice  that,  from  Proposition  3.4,  the  cumulants  of  satisfy  the  assumption 
of  Proposition  3.5.  Assume  that  cind  gi  are  the  polynomials  of  degrees  two  and 
three,  respectively,  in  Proposition  3.5.  Define  T„_i,  =  gi{tn,b)  and  T*;,  =  g2[tn,b)- 

Again,  from  Propositions  3.4  and  3.5,  one  can  see  that  /ci(t„,i,)  = 

«2(<n.t)  =  1  +  /C3(^t»,6)  =  0(m-*/2),  /C4(t„,6)  =  0(m-*/*);  Ki{Tn,b)  - 

0(m-^),  K2{Tn,b)  =  l+0(m"'),  Kz{Tn,b)  =  0(m"'),  K4{Tn,b)  =  Kl(r„*i,)  = 

0(m“M,  K2{T;i,)  =  1  -H  C>(m-^),  K3{T’^^)  =  0(m~^),  and  K4(r;,i,)  =  0(m"i).  Thus 
in  terms  of  the  distance  of  the  differences  between  the  first  four  cumulants  of  two 
random  variables,  both  Tn.i  and  T*  j,  converge  to  the  standard  normal  random  vari¬ 
able  fcister  than 

In  general,  we  need  to  estimate  the  cumulants  of  For  convenience,  we 
denote,  for  i  =  1,2, 3, 4, 

K,  =  K,(t„,6),  (3-14) 

a.s  an  estimate  for  the  ith  cumulant  of  Following  the  terminology  in  Glynn 
[12],  we  define 

tn.b  =  (X-EX)/{Vr,,k/ny/^  (3.15) 
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as  the  zeroth  order  pivot  (the  usual  approach); 


Tn,b  =  5l(<n,6)  =  +  («3/6  '  «l)  "  [{^^2  "  l)/2]t„.i  -  {K3/6)tl  i,  (3.16) 

as  the  first  order  pivot;  and 

^  =  52^6) 

=  tn,b  +  («^3/6  —  Ki)  +  [K4/8  —  {k2  —  l)/2]f„,i, 

-(«3/6)^^,6-(«4/24)<,  (3.17) 

as  the  second  order  pivot.  Notice  that  is  the  pivot  associated  with  the  traditional 
batch  means  method;  r„,6  is  an  estimate  of  which  is  the  unique  polynomial  of 
degree  two  of  in  the  sense  of  Proposition  3.5;  and,  similarly,  ^  is  an  estimate 
of  T’h,  which  is  the  corresponding  unique  polynonoial  of  degree  three  of  t.  Thus  if 

the  estimators  ^,’s  are  reasonably  well-behaved,  one  can  expect  that  both  Tn,t,  and 
r*  (,  converge  to  the  standard  normal  random  variable  faster  than 

To  construct  confidence  intervals  it  is  computationally  more  convenient  to  use 
the  inverted  Cornish-Fisher  expansions  (Hall  [12]).  For  the  ^-quantile  point,  zs,  of 
the  standard  normal  distribution  function,  the  100(1  —26)%  confidence  intervals 
for  the  three  pivots  for  the  batch  means  method  are 

[X,^  -  h'{zs){Vr.,blnY'\Xm  -  h'{-zs){Vr,,, 1X1^1%  (3.18) 

where  h'{z)  =  z  for  h'{z)  =  z  -h  (/Tj  -  /Ta/G)  -|-  [(/r2  -  l)/2]2  +  {k^I^)z'^  for 
and  h'{z)  =  z  -f  («i  —  ks/G)  +  [(«:2  ~  l)/2  —  («:4/8)]2:  -i-  (ks/G)*^  -|-  («:4/24)c^  for 

It  follows  from  equation  (3.18)  that  the  lengths  of  these  100(1  —  26)7c  con¬ 
fidence  intervals  are  as  follows:  2zi(  V'„.i,/n)*^^  for  the  zeroth  order  pivot  tny  2z((  1  -f 
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{k2  -  l)/2)(V'„,(,/n)^''^  for  the  first  order  pivot  Tny,  and  [(l+(<r2  -  l)/2-(«^/8))2zf+ 
(/r4/12)2i](V'„,i,/n)*/^  for  the  second  order  pivot  r*^.  We  have  shown  in  Proposi¬ 
tion  3.4  that  K2  —  1  =  0(n~^)  +  0{b~^)  and  K4  =  0{n~^).  Thus  if  both  estimators 
Kj  and  K4  are  reasonably  well-behaved,  then  one  can  expect  that  the  increase  of 
the  length  of  the  confidence  intervals  generated  by  the  Johnson-Glynn  pivots  are 
asymptotictdly  negligible. 

3.5.  ComputationaJ  Efficiency- 

In  this  section  we  will  compare  the  amount  of  computation  required  for  the 
traditional  batch  means  method,  first  and  second  order  Johnson-Glynn  pivots  for 
batch  means  method,  and  the  regenerative  method  of  simulation. 

\  direct  comparison  of  the  first  two  methods  with  the  regenerative  method  is 
difficult.  For  our  purpose,  we  assume  that  there  are  m  samples,  which  are  divided 
into  n  batches  of  b  samples  each,  and  there  are  uq  complete  regenerative  cycles 
within  these  m  samples.  Note  that,  by  strong  law  of  large  number,  tiq  w  mf  Et, 
where  Er  is  the  expected  length  of  the  regenerative  cycles. 

We  have  the  following  observations. 

(1)  Generating  m  samples:  for  many  simulation  studies  of  real  world  processes, 
the  amount  of  computation  required  of  generating  m  samples  is  of  the  form 
m  •  Co  +  0(1),  where  co  is  a  constant  which  depends  on  the  real  world  system 
been  simulated,  implementation  of  the  simulation  program,  and  the  computer 
system  used. 

(2)  Computing  batch  means:  the  amount  of  computation  needed  is  m  •  Ci  +  0(  1 ). 
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where  ci  does  not  dependent  on  m. 

(3)  Computing  sample  mean:  there  are  two  Ccises,  first,  if  the  sample  mccin  is 
computed  directly  from  the  samples  then  the  amount  of  computation  needed 
is  m  •  C2  +  0(1);  on  the  other  hand,  if  the  saunple  mean  is  computed  after 
the  batch  means  are  generated  then  the  amount  of  computation  needed  is 
n  ■  C2  +  0(1),  where  C2  does  not  depend  on  either  m  or  n. 

(4)  Computing  sample  variance:  the  «imount  of  computation  needed  is  m  •  ca  + 
0(1),  where  ca,  similar  to  cq,  does  not  depend  on  m. 

(5)  Computing  K4  for  Johnson-Glynn  pivots:  the  amount  of  computation 

needed  is  n -04 +  0(1),  n -05  + 0(1),  and  n -06  +  0(1),  respectively,  where  C4,  C5. 
and  ce  are  independent  of  n  but  are  dependent  on  the  system  been  simulated, 
implementation  of  the  simulation  program,  and  the  computer  system  used. 

(6)  Computing  V^’s  and  r.’s  for  regenerative  simulation:  the  amount  of  computa¬ 
tion  needed  is  m  •  C7  +  0(1),  and  m  •  cg  +  0(1),  respectively,  w'here  cr  and  cg 
are  constants. 

(7)  Computing  Xm  and  r  for  regenerative  simulation:  the  amount  of  computation 
needed  is  no  •  cg  +  0(1),  and  no  •  cio  +  0(1),  respectively,  where  cg  and  cio  are 
constants. 

(8)  Computing  for  regenerative  simulation:  the  amount  of  computation  needed 
is  no  •  Cji  +  0(1),  where  Cu  does  not  depend  on  no- 

(9)  It  is  easy  to  show  that  constants  ci  =  C2  =  c?  =  Cg  =  Cg  =  Ciq.  Assume 
that  the  operations  of  ith  power  of  a  constant  will  take  the  same  amount  of 
computational  effort  for  f  =  2,  3,  and  4,  then  we  have  cg  =  C4  =  C5  =  ce. 
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We  have  the  following  results. 


Proposition  3.7.  Under  the  conditions  and  constants  specified  in  Observations 

(1)  to  (9)  above,  the  amount  of  computation  required  for  confidence  interval  meth¬ 
ods  of  simulation  output  analysis  is  as  follows. 

(1 )  Traditional  i.i.d.  method:  m(cQ  +  Ci  +  cj)  +  0(1). 

(2J  Traditional  batch  means  method:  m(co  +  ci)  +  n(ci  +  C3)  -(-  0(1). 

(3)  First  order  Johnson-Glynn  pivot  for  the  batch  means  method:  m(co  +  cj)  + 
n(ci  +  2C3)  +  0(1). 

(4)  Second  order  Johnson-Glynn  pivot  for  the  batch  means  method:  m(co  +  ci)  + 
n(ci  +  3C3)  +  0(1). 

(5)  regenerative  method:  m(co  +  2ci)  +  no(2ci  -f  cn)  +  0(1). 

Proof:  (1)  Traditional  i.i.d.  method  needs  to  generate  samples,  and  compute  sam¬ 
ple  mean  and  sample  variance. 

(2)  Traditional  batch  means  method  needs  to  generate  samples,  and  compute  batch 
means,  sample  mean,  and  kj. 

(3)  First  order  Johnson-Glynn  pivot  for  the  batch  means  method  needs  to  generate 
samples,  and  compute  batch  means,  sample  mean,  K2,  and  K3. 

(4)  Second  order  Johnson-Glynn  pivot  for  the  batch  means  method  needs  to  gen¬ 
erate  samples,  and  compute  batch  means,  sample  mean,  K2,  «3,  and  K4. 

(5)  Regenerative  method  needs  to  generate  samples,  and  compute  i'^’s,  rjs,  Xm,  T, 

and  □ 

It  can  be  seen  that  the  amount  of  computation  required  for  the  Johnson-Glynn 
pivots  for  the  batch  means  method  is  more  than  that  of  the  traditional  batch  means 
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method.  However,  the  relative  increments  are  of  0(b  ‘),  which  is  asymptotically 
negligible  as  b  —*  cx),  for  either  cases. 
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Chapter  4 


Numerical  Results 


In  this  chapter  we  report  the  results  of  some  Monte  Carlo  studies  for  the 
coverage  statistics  of  “normal  quantile”  confidence  intervals  based  on  the  Johnson- 

Glynn  pivots  t„.i,  r„,6,  and  T’y 


4.1.  Notation  and  Precision 

Let  us  define  Ii,  1^,  h  as  follows:  Ii  represents  the  fraction  of  rephcations  for 
which  the  exact  value  Ues  to  the  left  of  the  confidence  interval;  I2  represents  the 
fraction  that  lies  in  the  confidence  interval;  and  I3  represents  the  fraction  that  lies 
to  the  right  of  the  confidence  interval.  Thus  I2  is  the  usual  coverage  fraction,  and 
Ii  and  I3  cure  the  one-sided  coverage  probabilities. 

For  each  of  our  examples,  we  make  2500  independent  replications  and  report 
empirical  coverage  fractions  Ii,  I2,  h,  sample  mean  of  the  length  of  the  confi¬ 
dence  interval  (SM),  seimple  standard  deviation  of  the  length  of  the  confidence 
interveil  (SSD),  and  sample  coefficient  of  variation  of  the  length  of  the  confidence 
intervad  (SCV),  which  is  the  ratio  of  sample  standard  deviation  over  sample  mean 
(SSD/SM). 

Notice  that  the  empirical  coverage  fractions  are  essentiailly  the  sample  means 
of  i.i.d.  binomial  random  variables  with  suitable  parameter  p.  For  2500  rephcations 
and  a  90%  confidence  interval,  li  and  I3  are  the  sample  means  of  2500  i.i.d.  binomial 
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random  variables  with  p  approximately  equal  to  0.05  (5%),  which  has  standard 
deviation  about  0.00436  (0.44%).  Similarly,  for  I2  the  p  value  of  the  corresponding 

1.1. d.  binomial  random  variables  is  approximately  0.90  (90%),  thus  the  standard 
deviation  of  the  stimple  mean  is  0.006  (0.6%).  These  are,  we  feel,  acceptable  levels 
of  accuracy  for  such  experiment. 

4.2.  Examples 

Example  4.1.  Let  be  i.i.d.  random  variables  with  distribution  func¬ 
tion  P{Yi  >  y)  =  1  —  y  <  0;  0,  otherwise.  (We  remark  parenthetically 

that  this  distribution  function  happens  to  be  the  stationary  distribution  of  the 
waiting  times  in  an  M/M/1  queue  with  mean  interarrival  time  1/A  and  mean  ser¬ 
vice  time  lip.)  See  Table  1  for  the  empirical  results  of  the  coverage  rate  for  five 
different  methods  of  output  analysis. 

Example  4.2.  Let  be  i.i.d.  random  variables  with  distribution  func¬ 
tion  P(Yi  >  y)  =  y  >  —1;  0,  otherwise.  Thus  each  Y  is  a  centered- 

exponential  random  variable  with  parameter  1.  This  example  was  studied  in  Efron 
[11]  and  Glynn  [12].  See  Table  2  for  the  empirical  results. 

Example  4.3.  Let  Xi,. be  samples  from  an  autoregressive  model  of  order 
one  such  that,  for  each  i,  Afi+i  =  0.5A',-  +  €,+1,  where  the  residuals  Ci ’s  are  centered- 
exponential  random  variables  discussed  at  Example  4.2.  This  example  was  studied 
in  Titus  [25].  See  Tables  3  and  4  for  the  empirical  results. 

Example  4.4.  Let  {Wi  :  t  >  1}  6e  a  sequence  of  waiting  times  in  an  M/M/1 
queue  with  arrival  rate  A  =  0.5  and  service  rate  p  =  1.  This  example  was  also 
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studied  in  Glynn  [12]  and  Titus  [25].  To  see  the  eiFect  of  initial  conditions,  we 
choose  several  different  initial  random  variable  Wi.  Let  W  be  the  random  variable 
distributed  as  the  stationary  waiting  times.  Those  cases  we  consider  include:  Wi  = 
0;  Wi  =  W ;  Wi  =  EW;  and  Wi  =  2EW.  See  Tables  5  to  12  for  the  empirical 
results. 


4.3.  Discussion  of  Numerical  Results 

Note  that  in  our  numerical  examples,  essentially  in  every  case  there  is  some 
improvement  in  actual  coverage  fraction  I2  from  zeroth  order  pivot  (the  traditional 
batch  means  method)  to  the  first  and  second  order  pivots.  As  a  matter  of  fact, 
Table  1  shows  that  for  our  Example  1,  the  first  order  Johnson-Glynn  pivot  com¬ 
petes  favorably  with  the  method  of  known  variance.  In  addition,  the  first  and 
second  order  pivots  tend  to  balance  the  one-sided  coverage  probabilities  /j  and  I3, 
moving  them  towards  their  desired  values  of  0.05  (5%).  This  confirms  the  asym¬ 
metry  corrections  induced  by  the  Johnson-Glynn  pivotal  transformations.  On  the 
other  hand,  we  note  that  a  confidence  interval  with  balanced  one-sided  coverage 
probabilities  does  not  necessarily  have  the  shortest  possible  length. 

One  may  notice  that  in  our  simulation  results,  the  second  order  Johnson- 
Glynn  pivot  only  provides  a  level  of  coverage  fraction  comparable  with  that  of  the 
first  order  pivot.  One  reason  for  this  is  that  in  our  analysis,  we  use  the  formal 
Cornish-Fisher  expansion  for  1^,6  in  deriving  both  pivots  T„_4  and  In  this 

way,  T*  j  is  a  function  of  in,b  and  its  sample  cumulants,  instead  of  that  of  the  first 
order  pivot  Tn^b-  Hence  the  T*  j,  we  derive  has  a  potential  of  improvement  over 
instead  of  that  of  T„,6.  There  is,  however,  another  approach  in  deriving  the  second 
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order  pivot.  If  we  can  derive  the  formal  Cornish-Fisher  expansion  of  J’n.j,  then  the 
second  order  pivot  T*  ^  is  then  a  function  of  Tnj,  and  its  sample  cumulants.  By  the 
same  idea  the  first  order  pivot  Tn,;,  corrects  some  effects  of  estimation  errors  and 
the  skewness  effect  associated  with  may  correct  those  estimation  errors 

and  the  skewness  effect  in  the  estimation  of  the  first  order  pivot.  Although  the 
second  approach  is  preferable  from  a  theoretical  point-of-view,  the  derivation  of  an 
additional  pivot  from  r„,6  will  be  much  more  difficult. 

To  extimine  the  effects  of  nonstationarity  and  the  impact  of  different  initial 
conditions,  we  choose  several  initial  random  variable  W\  in  the  simulation  of  Exeun- 
ple  4.4.  Let  W  be  the  random  v'ariable  distributed  as  the  stationary  waiting  times. 
Those  cases  we  consider  include:  Wy  =  0;  ILi  —  W;  Wx  =  EW\  and  VFi  =  2E\V . 
As  we  can  see  empirically  from  Tables  5  to  12,  for  the  waiting  times  of  an  M/M/1 
queue,  both  the  nonstationarity  and  the  different  initial  conditions  have  only  a  very 
insignificant  etfect  on  the  coverage  fraction,  one-sided  coverage  probabilities,  and 
the  length  of  the  confidence  intervals.  Further  efforts  will  be  needed  to  verify  this 
theoretically. 

We  can  also  observe  that  the  Johnson-Glynn  pivots  produce  longer  confidence 
intervals  on  average,  and  these  confidence  intervals,  in  general,  are  more  variable 
than  those  of  the  traditional  batch  means  method.  However,  as  shown  in  Chapter 
3,  the  increase  of  length  in  confidence  intervals  is  asjTnptotically  negligible  as  n 
and  h  increase.  Moreover,  due  to  the  fact  that  many  confidence  intervals  produced 
from  the  traditional  batch  means  method  do  have  an  undercoverage  problem,  this 
increase  of  length  seems  to  be  a  necessity  rather  than  a  liability.  The  more  vari¬ 
able  confidence  intervals  produced  by  the  Johnson-Glynn  pivots  are  also  a  natuial 
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property  associated  with  these  correction  methods.  Since  the  correction  terms  cire 
stochastic,  in  general  more  constcints  need  to  be  estimated  and  then  additional 
vcuiance  is  introduced. 
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Table  1.  Coverage  Fraction  for  Example  4.1. 


normal 

t 

first 

second 

known  variance 

5 

0.6604 

0.7620 

0.7740 

0.7904 

0.9312 

10 

0.7588 

0.7912 

0.8512 

0.8472 

0.9180 

15 

0.8092 

0.8260 

0.8928 

0.8892 

0.9152 

20 

0.8324 

0.8464 

0.8992 

0.8956 

0.9080 

25 

0.8412 

0.8532 

0.9012 

0.8976 

0.8988 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  known  variance:  batch  means  method  with  known  variance 

•  6=1,  since  independently  identical  observations 

•  Coverage  fraction  for  the  first  order  pivot  competes  favorably  with  that  of  the 
batch  means  method  with  known  variance 
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Table  2.  Coverage  Fraction  for  Centered-Exponential  Random  Variables. 


Sample 

Coverage  Fraction 

Length  of  C.I. 

Size  (m) 

Pivot 

h 

h 

h 

SM 

SSD 

SCV 

normal 

1^^  1  in 

0.7140 

1.12 

Bra 

0.58 

t 

b|R  I 

0.7864 

IPs 

1.45 

ill 

0.58 

5 

first 

0.8220 

1.61 

1.01 

0.63 

second 

0.0160 

0.8188 

0.1652 

1.58 

0.98 

0.62 

normal 

0.7904 

0.1788 

BEOI 

0.39 

0.42 

t 

0.8252 

0.1560 

Mm 

0.43 

0.42 

10 

first 

lS][Syii9 

0.8660 

0.1124 

Ira 

0.58 

0.49 

second 

0.0216 

0.8652 

Hal 

0.56 

0.48 

normal 

0.8276 

0.1448 

Bra 

0.34 

t 

0.8488 

0.1336 

ira 

ira 

0.34 

15 

first 

0.8884 

0.0852 

in 

0.41 

second 

0.8880 

0.0856 

0.94 

0.38 

0.40 

normal 

0.8548 

0.1212 

0.69 

Bra 

Bra 

t 

0.8668 

0.1140 

0.73 

ira 

ira 

20 

first 

0.8992 

0.0764 

0.82 

ira 

second 

0.0264 

0.8952 

0.0784 

0.81 

0.28 

ira 

normal 

0.8632 

0.1144 

0.63 

Bra 

Bra 

t 

■IKtIKf 

0.8744 

0.1072 

0.65 

ira 

ira 

25 

first 

0.0272 

0.8984 

0.0744 

0.72 

ira 

second 

0.0284 

0.8968 

0.0748 

0.72 

ill 

0.32 

•  2500  independent  replications 

•  SM:  sample  mean 

•  SSD:  sample  standard  deviation 

•  SCV:  sample  coefficient  of  variation 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction 

•  One-sided  coverage  probabilities  /j  and  Iz  of  the  first  and  second  order  pivots 
are  more  balanced 
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Table  3.  Coverage  Fraction  for  An  AR(1)  Process. 


Coverage  Fraction 

Length  of  C  .I. 

Pivot 

h 

I2 

h 

SM 

SSD 

SCV 

normal 

rfiwi  vft  M 

0.7900 

0.1312 

0.36 

0.055 

0.16 

t 

0.7936 

0.1296 

0.36 

0.056 

0.16 

2 

first 

0.8000 

0.1068 

0.37 

0.063 

0.17 

second 

0.7992 

0.1072 

0.37 

0.062 

0.17 

normal 

0.0568 

0.8364 

HjY 

0.34 

0.066 

BSSi 

i 

0.0496 

0.8524 

■on  ? :  9 

0.35 

0.068 

im 

4 

first 

0.0608 

0.8592 

Bli :  i  9 

0.37 

0.084 

im 

second 

0.0620 

0.8580 

EKs!9 

0.37 

0.082 

normal 

0.8460 

teiiMiJI 

t 

Kinlf39 

0.8664 

Int  >  ‘  9 

im 

Bfiwa 

6 

first 

0.0488 

0.8764 

0.0748 

0.36 

BE^a 

0.25 

second 

0.0520 

0.8728 

0.0752 

0.36 

0.089 

0.25 

normal 

0.8624 

ten 

fe'VlM 

t 

BISk  ^I 

0.8788 

iis 

BE^I 

im 

8 

first 

Bm 

0.8912 

BKi  (^9 

0.36 

0.10 

19 

second 

BfifEW 

0.8892 

0.0748 

0.36 

0.099 

0.28 

normal 

0.8492 

knutii  9 

WEM 

19 

t 

0.8832 

Bn^9 

BfilJB 

19 

10 

first 

BtiRaM 

0.8968 

0.0652 

iia 

BnuB 

0.30 

second 

0.8948 

0.0664 

0.35 

0.10 

0.29 

0.0568 

0.8312 

■n  IF  9 

BliSB 

19 

■nil 

0.0356 

0.8812 

Bifi?9 

i9 

15 

0.0340 

0.8944 

0.36 

0.12 

im 

second 

0.0344 

0.8920 

0.0680 

0.35 

0.11 

19 

•  A,+i  =  0.5  Xi  +  e^+i 

•  e,:  centered-exp(l) 

•  =  0 

e  5’mple  size:  m  =  120 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 


•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction 

•  One-sided  coverage  probabilities  Ij  and  I3  of  the  first  and  second  order  pivots 
are  more  bcilanced 
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Table  4.  Coverage  Fraction  for  An  AR(1)  Process. 


Batch 

Coverage  Fraction 

Length  of  C.I. 

Size  (b) 

Pivot 

h 

h 

h 

SM 

SSD 

SCV 

normal 

0.0860 

0.7884 

0.1256 

0.26 

0.029 

0.11 

t 

0.0S56 

0.7904 

0.1240 

0.26 

0.029 

0.11 

2 

first 

0.7924 

0.1076 

0.26 

0.032 

0.12 

second 

0.1004 

0.7920 

0.1076 

0.26 

0.031 

0.12 

normal 

0.8496 

0.0936 

0.25 

0.034 

0.14 

t 

Bnlf  ( « 

0.8548 

0.0916 

0.25 

0.034 

0.14 

4 

first 

0.0644 

0.8580 

0.0776 

0.26 

0.039 

0.15 

second 

0.0648 

0.8572 

0.0780 

0.26 

0.038 

0.15 

normal 

0.8628 

0.0884 

0.24 

0.037 

0.16 

t 

■infill 

0.8696 

0.0860 

0.25 

0.038 

0.16 

6 

first 

■tnlilM 

0.8688 

0.0764 

0.25 

0.045 

0.18 

second 

0.0556 

0.8676 

0.0768 

0.25 

0.044 

0.17 

normal 

0.8664 

0.0848 

0.23 

0.040 

0.17 

t 

0.8768 

0.0788 

0.24 

0.042 

0.17 

8 

first 

■inifvl 

0.8832 

0.0692 

0.25 

0.050 

0.20 

second 

0.0492 

0.8816 

0.0692 

0.25 

0.049 

0.20 

normal 

0.8668 

0.0864 

0.23 

0.041 

0.18 

t 

anvils 

0.8828 

0.0784 

0.24 

0.043 

0.18 

10 

first 

Hjffi '  Ifl 

0.8904 

0.0656 

0.25 

0.051 

0.21 

second 

■wit  1 3 

0.8892 

0.0656 

0.25 

0.050 

0.20 

normal 

vnb  ifiii 

0.8652 

0.0812 

0.22 

0.047 

0.22 

t 

nntn 

0.8856 

0.0732 

0.23 

0.050 

0.22 

15 

first 

0.0452 

0.8884 

0.0664 

0.24 

0.060 

0.24 

second 

0.0456 

0.8876 

0.0668 

0.24 

0.058 

0.24 

•  -Yx+i  —  0.5  A,  + 

•  e,:  centered-exp(l) 

•  Ai  =  0 

•  sample  size:  m  =  240 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction 

•  One-sided  coverage  probabilities  I\  and  of  the  first  eind  second  order  pivots 
are  more  baleinced 
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Table  5.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Coverage  Fraction 


Batch 
Size  (b) 

Pivo: 

20 

normal 

t 

first 

second 

25 

normal 

t 

first 

second 

40 

normal 

t 

first 

second 

50 

normal 

t 

first 

second 

100 

normal 

i 

first 

second 

125 

normal 

t 

first 

second 

0.0608 

0.0644 


0.0520 

0.0540 


0.0356 

0.0292 

0.0400 

0.0416 


0.0440 


0.0388 

0.0260 

0.0284 

0.0296 


0.8104 

0.8200 

0.8408 

0.8368 


0.8208 

0.8268 

0.8568 

0.8524 


0.8272 

0.8400 

0.8668 

0.8648 


0.8212 

0.8408 

0.8608 

0.8584 


0.8124 

0.8492 

0.8764 

0.8744 


0.7980 

0.8432 

0.8672 

0.8636 


536 
268 
0.1016 
0.1032 


•  M/M/1  waiting  limes 

•  /?  =  0.5 

•  U'l  =  \V 

•  {1^/  :  i  >  1}  is  stationary 

•  sample  size:  m  =  1000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  /n 

•  One-sided  coverage  probabilities  Ii  and  I3  of  the  first  and  second  order  pivots 
are  more  balajiced 
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Table  6.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Batch 

Coverage  Fraction 

Length  of  C.I. 

Size  (6) 

Pivot 

h 

h 

h 

SM 

SSD 

scv 

normal 

■ITOin^ 

0.8224 

|TO  inn 

0.32 

0.082 

0.25 

i 

■SVSfMfl 

0.8256 

■nt 

0.33 

0.083 

0.25 

20 

first 

0.8372 

0.36 

0.10 

0.29 

second 

0.8348 

■nfi'  tifeB 

0.35 

0.10 

0.29 

normal 

0.8336 

0.27 

t 

■SKm 

0.8348 

IsImS 

■iKgM 

VatI'I 

0.27 

25 

first 

0.8512 

0.0920 

0.37 

■OlrH 

0.32 

second 

0.0580 

0.8500 

0.0920 

0.37 

0.11 

0.31 

normal 

0.35 

0.30 

t 

0.36 

0.30 

40 

first 

0.0492 

0.8652 

0.39 

0.38 

second 

0.0504 

0.8624 

0.0872 

0.39 

■IKtiM 

normal 

0.0300 

0.8484 

0.1216 

Bia 

Em 

0.31 

t 

0.0292 

0.8516 

0.1192 

loiB 

0.31 

50 

first 

0.0448 

0.8692 

0.0860 

0.40 

0.39 

second 

0.0468 

0.8660 

0.0872 

0.40 

lilU 

0.38 

normal 

0.8480 

0.36 

Em 

0.34 

t 

0.8620 

0.38 

0.34 

100 

first 

0.0340 

0.8788 

0.0872 

0.42 

0.18 

0.42 

second 

0.0380 

0.8748 

0.0872 

0.41 

0.17 

0.41 

normal 

0.8328 

0.1296 

0.35 

0.36 

i 

0.8516 

0.1196 

0.38 

EHfl 

0.36 

125 

first 

0.8720 

0.0924 

0.42 

0.19 

0.46 

second 

0.0360 

0.8712 

0.0928 

0.42 

0.19 

0.44 

•  M/M/1  waiting  times 

•  p  =  0.5 

9  \\\  =  W 

•  {W,  :  I  >  1}  is  stationary 

•  sample  size;  rn  =  2000 

•  2500  independent  replications 

•  normal;  batch  means  method  with  normal  quantiles 

•  t;  batch  means  method  with  t  quantiles 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  Ii  and  I3  of  the  first  and  second  order  pivots 
are  more  balanced 
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Table  7.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Coverage  Fraction 

Length  of  C.I. 

Eiil 

Pivot 

h 

I2 

h 

SM 

SSD 

SCV 

normal 

0.8068 

0.45 

Bsa 

t 

Bnn  ( ! » 

0.8152 

0.46 

Bm 

20 

first 

0.8428 

■iMlIiK  !■ 

0.51 

ngfl 

0.40 

second 

0.0624 

0.8380 

0.0996 

0.51 

0.39 

normal 

0.8192 

ncfl 

0.35 

t 

0.8268 

iSi 

nM 

0.35 

25 

first 

■W  iflllfl 

0.8532 

!|B« 

0.43 

second 

0.0544 

0.8480 

BH 

nEM 

0.42 

normal 

0.0336 

0.8216 

Esa 

2a 

0.39 

t 

0.0284 

0.8336 

olll 

0.39 

40 

first 

0.0392 

0.8648 

0.0960 

Mm 

0.50 

second 

0.0412 

0.8620 

0.0968 

0.54 

0.49 

normal 

■iliism 

0.8204 

0.1460 

0.48 

SSI 

Bsa 

t 

0.8344 

0.1368 

0.50 

ilafl 

50 

first 

0.0392 

0.8604 

0.1004 

0.58 

Sl@l 

HI 

second 

0.0412 

0.8580 

0.1008 

0.57 

0.28 

BH 

normal 

EBI 

WEM 

i 

Bm  1 

■IIWIgM 

HI 

100 

first 

0.8752 

0.0968 

BH 

B^ 

BH 

second 

0.0284 

0.8732 

0.0984 

0.59 

0.31 

0.52 

normal 

0.0440 

0.7968 

■iKiril 

0.22 

Bsa 

t 

0.0272 

0.8412 

■aMO 

0.25 

125 

first 

0.0296 

0.8632 

0.1072 

0.60 

0.34 

IH 

second 

0.0316 

0.8592 

0.1092 

0.60 

0.33 

BH 

•  M/M/1  waiting  times 

•  p  =  Q.b 

•  \\\  =  0 

•  :  i  >  1}  is  only  asymptotically  stationary 

•  sample  size:  m  =  1000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second;  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  1^ 

•  One-sided  coverage  probabilities  l\  and  I3  of  the  first  and  second  order  pivots 
are  more  balanced 
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Table  8.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Batch 

Coverage  Fraction 

Length  of  C.I. 

Size  (5) 

Pivot 

h 

I2 

h 

SM 

SSD 

scv 

normal 

0.8216 

■iTiw9 

■3^9 

t 

0.8228 

kk3 

BE  9 

20 

first 

0.0640 

0.8372 

BX$>9 

0.36 

0.10 

0.29 

second 

0.0660 

0.8340 

0.1000 

0.35 

0.10 

0.29 

normal 

0.8336 

WEM 

im 

t 

■ijfl  m  M 

0.8356 

Wm 

25 

first 

Bijfl  til 

0.8520 

0.0920 

0.37 

im 

second 

0.0572 

0.8500 

0.0928 

0.37 

am 

ifa 

normal 

0.0332 

0.8416 

0.1252 

0.35 

0.30 

t 

0.0308 

0.8472 

0.1220 

0.36 

0.30 

40 

first 

0.0464 

0.8672 

0.0864 

0.39 

0.38 

second 

0.0476 

0.8652 

0.0872 

0.39 

EiEa 

0.37 

normal 

ii*h*p|*  w 

0.8456 

0.11 

t 

0.8516 

0.1212 

na 

0.11 

bsi 

50 

first 

0.8680 

0.40 

0.16 

second 

0.0452 

0.8840 

0.0884 

0.40 

0.15 

Bm 

normal 

0.8472 

0.1236 

0.36 

KUH 

0.34 

t 

0.8628 

0.1144 

0.37 

0.34 

100 

first 

Bju?  1 9 

0.8760 

0.42 

0.18 

0.42 

second 

0.0352 

0.8724 

0.0924 

0.41 

0.17 

0.41 

normal 

0.8324 

0.35 

EiEa 

0.36 

t 

■Sill'S  9 

0.8520 

0.38 

0.36 

125 

first 

0.0336 

0.8704 

0.0960 

0.42 

0.19 

0.46 

second 

0.0336 

0.8700 

0.0964 

0.42 

0.18 

0.44 

•  M/M/1  waiting  times 

•  p  =  0.5 

•  =  0 

•  {W,  :  z  >  1}  is  only  asymptotically  stationary 

•  sample  size:  m  =  2000 

•  2500  independent  replications 

•  normal:  batch  meeins  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  /j  and  Iz  of  the  first  and  second  order  pivots 
are  more  balanced 


Table  9.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Coverage  Fraction 


Pivot 


normal 

t 

first 

second 


normal 

t 

first 

second 


normal 

t 

first 

second 


normal 

t 

first 

second 


normal 

t 

first 

second 


normal 

t 

first 

second 


Length  of  C.I. 


SM  SSD  SCV 


0.0440 

0.0400 

0.0600 

0.0636 


0.0404 

0.0428 


0. 

0. 

0. 
0.0296 


0.0456 

0.0280 

0.0312 

0.0324 


0.8092 

0.3212 

0.8440 

0.8392 


0.8220 

0.8288 

0.8516 

0.8472 


0.8256 

0.8348 

0.8672 

0.8640 


0.8232 

0.8388 

0.8608 

0.8576 


0.8108 

0.8492 

0.8780 

0.8748 


0.7972 

0.8436 

0.8652 

0.8628 


.15 

0.33 

.15 

0.33 

.20 

0.40 

.20 

0.39 

.16 

0.35 

.17 

0.35 

.23 

0.42 

.22 

0.41 

.19 

0.39 

.19 

0.39 

.28 

0.50 

.27 

0.49 

0. 

.0. 

0. 
0.0956 


572 
284 
036 
0.1048 


•  M/M/1  waiting  times 

•  ^  =  0.5 

•  =  EW 

•  {Wi  :  i  >  1}  is  only  asymptotically  stationary 

•  sample  size;  m  =  1000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson- Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  Ii  and  I3  of  the  first  and  second  order  piv'ots 
are  more  balanced 
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Table  10.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Batch 

Coverage  Fraction 

Length  of  C.I. 

Size  (6) 

Pivot 

h 

h 

SM 

SSD 

scv 

normal 

0.8212 

0.1320 

0.33 

13^ 

t 

0.8240 

0.1308 

0.33 

20 

first 

0.8388 

0.0964 

0.36 

0.10 

0.29 

second 

0.0664 

0.8364 

0.0972 

0.35 

0.10 

0.29 

normal 

0.8332 

■t  w.  i  w 

wnm 

BIS 

t 

0.8352 

mmiM 

25 

first 

0.0572 

0.8516 

0.37 

■SVrH 

0.32 

second 

0.0580 

0.8500 

0.0920 

0.37 

0.11 

0.31 

normal 

0.8436 

0.30 

t 

■WKitf 

0.8488 

ims 

ii9 

0.30 

40 

first 

0.0476 

0.8668 

0.0856 

0.39 

0.38 

second 

0.0484 

0.8652 

0.0864 

0.39 

0.37 

normal 

0.8456 

0.1236 

0.35 

0.31 

t 

0.8524 

0.1200 

0.36 

0.31 

50 

first 

0.8712 

0.0860 

0.40 

0.39 

second 

0.8656 

0.0880 

0.40 

3^1 

0.38 

normaJ 

Miliii  M 

0.8476 

0.1220 

0.36 

0.12 

t 

0.8644 

0.1128 

0.37 

0.13 

Utl 

100 

first 

0.8760 

wliiMa 

0.42 

0.18 

Mm 

second 

0.0360 

0.8728 

0.41 

0.17 

0.41 

normal 

0.0352 

0.8332 

0.1316 

0.35 

QEB 

Bsa 

t 

0.8524 

0.1204 

0.38 

kIH 

125 

first 

0.0336 

0.8720 

0.0944 

0.42 

0.19 

0.46 

second 

0.0340 

0.8708 

0.0952 

0.42 

0.18 

0.44 

•  M/M/1  waiting  times 

•  =  0.5 

•  =  EW 

•  {W,  :  i  >  1}  is  only  asymptotically  stationary 

•  sample  size:  m  =  2000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  Ii  and  I3  of  the  first  and  second  order  pivots 
are  more  balanced 


Table  11.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Coverage  Fraction 

Length  of  C.I. 

Pivot 

h 

h 

h 

SM 

SSD 

SCV 

normal 

0.0452 

0.8164 

0.1384 

0.45 

0.15 

0.33 

t 

0.0412 

0.8248 

0.1340 

0.46 

0.15 

0.33 

20 

first 

0.0640 

0.8444 

0.0920 

0.52 

0.20 

0.40 

second 

0.0672 

0.8400 

0.0916 

0.51 

0.20 

0.39 

normal 

0.0392 

0.8244 

0.1364 

0.46 

0.16 

0.35 

t 

0.0356 

0.8316 

0.1328 

0.48 

0.17 

0.35 

25 

first 

0.0540 

0.8544 

0.0916 

0.54 

0.23 

0.42 

second 

0.0572 

0.8492 

0.0936 

0.53 

0.22 

0.41 

normal 

0.0356 

0.8276 

0.1368 

0.48 

0.19 

0.39 

t 

0.0308 

0.8396 

0.1296 

0.50 

0.19 

0.39 

40 

first 

0.0420 

0.8672 

0.0908 

0.57 

0.29 

0.50 

second 

0.0440 

0.8644 

0.0916 

0.56 

0.27 

0.49 

normal 

0.0384 

0.8240 

0.1376 

0.48 

0.19 

0.40 

t 

0.0300 

0.8424 

0.1276 

0.50 

0.20 

0.40 

50 

first 

0.0428 

0.8640 

0.0932 

0.58 

0.29 

0.50 

second 

0.0440 

0.8612 

0.0948 

0.57 

0.28 

0.49 

normal 

0.0400 

0.8152 

0.1448 

0.47 

0.21 

0.44 

t 

0.0260 

0.8536 

0.1204 

0.53 

0.23 

0.44 

100 

first 

0.0292 

0.8824 

0.0884 

0.60 

0.32 

0.53 

second 

0.0304 

0.8800 

0.0896 

0.59 

0.31 

0.52 

normal 

0.0476 

0.7996 

0.1528 

0.47 

0.22 

0.47 

t 

0.0308 

0.8468 

0.1224 

0.54 

0.25 

0.47 

125 

first 

0.0324 

0.8688 

0.0988 

0.61 

0.34 

0.56 

second 

0.0324 

0.8672 

0.1004 

0.60 

0.33 

0.55 

•  M/M/1  waiting  times 

•  p  =  0.5 

•  =  2  E\V 

•  {W,  :  i  >  1}  is  only  asymptotically  stationary 

•  sample  size;  m  =  1000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  and  I3  of  the  first  and  second  order  pivots 
are  more  balanced 
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Table  12.  Coverage  Fraction  for  M/M/1  Waiting  Times. 


Kssai 

Coverage  Fraction 

Length  of  C.I. 

Pivot 

I2 

I3 

SM 

SSD 

SCV 

normal 

0.8204 

0.1308 

0.33 

0.082 

0.25 

t 

0.8240 

0.1288 

0.33 

0.083 

0.25 

20 

first 

0.8396 

0.0940 

0.36 

0.10 

0.29 

second 

■iIiWiIm 

0.8348 

0.35 

0.10 

0.29 

normal 

0.8344 

BBI 

t 

MR 

0.8388 

MR  .*9 

WSm 

25 

first 

InR  iitiil 

0.8540 

MRkS9 

0.37 

0.12 

second 

0.0592 

0.8512 

0.0896 

0.37 

0.11 

iiy 

normal 

0.8436 

0.35 

0.30 

t 

loRtEn 

0.8496 

MQ^9 

0.36 

0.30 

40 

first 

0.0484 

0.8684 

0.39 

0.38 

second 

0.0496 

0.8668 

0.39 

0.37 

normal 

J 

0.8476 

0.31 

t 

MR  V  9 

0.8568 

MQ^9 

na 

0.31 

50 

first 

MR  [f  9 

0.8712 

MRctem 

0.40 

0.39 

second 

0.0472 

■IKftWf 

0.0856 

0.40 

0.38 

normal 

Kin;  n 

0.8500 

0.36 

0.12 

0.34 

t 

0.8664 

■Win'  9 

0.37 

0.13 

0.34 

100 

first 

MRn9 

0.8784 

0.42 

0.18 

0.42 

second 

0.0368 

0.8752 

Mfi^9 

0.41 

0.17 

0.41 

normal 

0.8348 

HsgUI 

0.35 

ElEM 

Esa 

t 

0.8564 

0.38 

125 

first 

■nRSeKM 

0.8732 

0.0920 

0.42 

0.19 

0.46 

second 

0.0356 

0.8724 

0.0920 

0.42 

0.18 

0.44 

•  M/M/1  waiting  times 

•  p  =  0.5 

•  W\=  2  EW 

•  {IVj  :  i  >  1}  is  only  asymptotically  stationary 

•  sample  size:  m  =  2000 

•  2500  independent  replications 

•  normal:  batch  means  method  with  normal  quantiles 

•  t:  batch  means  method  with  t  quantiles 

•  first:  first  order  Johnson-Glynn  pivot 

•  second:  second  order  Johnson-Glynn  pivot 

•  Both  first  and  second  order  pivots  improve  coverage  fraction  I2 

•  One-sided  coverage  probabilities  /j  and  I3  of  the  first  and  second  order  pivots 
are  more  balanced 
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Chapter  5 


Conclusion 


As  the  use  of  simulation  becomes  more  popular,  the  need  for  a  method  of 
simulation  output  analysis  that  can  be  applied  to  a  large  class  of  stochaistic  processes 
becomes  more  important. 

In  this  dissertation,  we  have  derived  a  uniqueness  property  of  Cornish-Fisher 
expansions,  the  formal  Edgeworth  expansion  for  the  batch  means  method,  Johnson- 
Glynn  pivots  and  the  associated  confidence  intervals  for  the  batch  means  method, 
and  detailed  procedures  of  implementation. 

Johnson-Glynn  pivotal  transformations  have  provided  a  new  wa>  of  generat¬ 
ing  confidence  intervals.  In  applying  this  approach  to  the  batch  means  method, 
they  appear  to  behave  well  empirically  and  seem  to  be  a  robust  procedure  for  the 
examples  in  Chapter  4.  To  verify  the  general  applicability  of  this  method,  many 
more  examples  should  be  run. 

However,  there  are  three  possible  drawbacks  from  a  practical  point-of-view: 
more  computation  time  needed,  longer  confidence  intervals  on  average,  and  more 
variable  intervals. 

On  the  other  hand,  as  shown  in  Chapter  3,  the  increase  of  computing  time  is 
relatively  small,  and  the  increase  of  length  in  confidence  interval  is  aisjTTiptotically 
negligible  as  n  and  b  increase.  Moreover,  due  to  the  fact  that  many  confidence 
intervals  do  have  an  undercoverage  problem,  this  increase  of  length  seems  to  be  a 
necessity  rather  than  a  habihty. 
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We  would  like  to  roint  out  some  potential  areas  for  future  research  and  de¬ 
velopment.  First,  although  the  assumptions  we  make  are  reasonable,  we  may  want 
to  relax  the  assumptions  we  made  on  stationairity,  either  by  including  the  initial 
condition  or  considering  eisymptotically  stationary  processes.  Second,  it  is  possible 
that  the  Johnson-Glynn  pivots  can  be  applied  to  overlapping  batch  means  method. 
Another  interesting  possibility  is  to  apply  the  same  procedure  to  the  ratio  esti¬ 
mator  of  a  weakly  regenerative  process.  Finally,  we  may  want  to  run  many  more 
numerical  examples  to  verify  the  general  apphcability  of  this  method. 
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Appendix 


Proof  of  Proposition  2.1;  (1)  Suppose  5i(X„)  =  +  a„,o  +  an,iXn  +  an,2A'^, 

where  a,  =  i  =  0,  1,  and  2.  Then  Ki((/i(X„))  =  EXn  +  an,o  +  o.n,\EXn  + 

^n^EX^  =  In.l  +  Qn.O  +  an.2  +  0(n~^),  K2(gi(X„))  =  1  +  2a„.i  +  /„,2  +  0(n~^), 

K3{gi(X„))  =  /„,3  +  6a„,2  +  0(n“^).  Notice  that  ki(^)  =0,  /C2(0  =  1,  and  K3(^)  =  0. 
If  we  need  \K^{()  —  «:,(5i(,Y„))|  =  0(n~^)  for  each  i,  1  <  f  <  3,  the  unique  solution 
of  0,’s  in  5i(A'’„),  up  to  order  is  such  that  a„,o  =  —U,\  +  ^n,3/6  + 

<2n,i  =  0{n~'^),  and  a„.2  =  -/n.a/d  +  0{n~^). 

(2)  Suppose  g2{Xn)  =  +  a„,o  +  a„,iX„  +  a„,2X^  +  a„,3X^,  and  a,  =  0{n~'^^^)  for 

t  =  0,  1,  2,  and  3.  Then  Ki(52(-Yn))  =  /n.i  +  o„.o  +  +  a„,2  +  a„,3/„,3  +  o(n"^), 

'^2(5'2(-^n))  =  1  +  2a„,i  T  +  /„,2  +  2a^_2  +  +  20n,2^n,3  +  4a„.2/„,i  +  6a„3  + 

6an,inn,3  +  n(n~^),  «:3(52(A^n))  =  ln,3  +  3a„,i/„;}  +  6a„,2  +  12a„,:a„,2  +  27a„,3/n,3  + 
lSnn.3^n,l  +  o(n~^),  «:4(52(-Vn))  =  ^n,4  +  24an,2ln,3  +  24an.3  +  72a„,ian,3  +  48a^2  + 
432a^3  +  o(n“').  If  we  need  iKi(^)  —  «t(52(A„))|  =  o{n~^)  for  each  i,  1  <  i  <  4, 
the  unique  solution  of  a,’s  in  g2{Xn),  up  to  order  (9(n~^)  is  such  that  Ono  = 

~L,l  +  ^n,3/6  4-  o(n~^),  a„,i  =  —ln,2l2  +  /n, 14,3/3  +  4,4/8  —  (7/36)42  + 

ar,,2  =  -4,3/6  +  o(n"^),  and  a„.3  =  -4,4/24  +  4,3/9  +  o{n~^).  □ 

Proof  of  Proposition  2.2:  (1)  Suppose  Ai(0  =  ^  +  n„,o  +  an,i<f  +  n„,2<f^,  where 
a,  =  i  =  0,  1,  and  2.  Then  ki(/ii(0)  =  an,o  +  a„,2,  K2(*i(0)  =  1  + 

2an.i  +  0(n"M.  K3{hi(0)  =  6a„,2  +  0(71-^).  Notice  that  ki(A'„)  =  0,  K2(A4)  =  1. 
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and  K3(.Y„)  =  /„,3  =  0(n~^/^).  If  we  need  |«i(A’„)  -  k,(/ii(0)I  = 

1,  1  <  i  <  3,  the  unique  solution  of  a,’s  in  hi{^),  up  to  order  0(n~^^^)  is  such  that 

fln.o  =  ln,i  -  ln,z/^  +  ^n,!  =  and  a„,2  =  /„.3/6  +  0(n~^). 

(2)  Suppose  h2{0  =  ^  +  an.o  +  «n,i^  +  On,2^*  +  o„,3^^,  and  a,  =  0(n"^/^)  for  f  =  0,  1, 

2,  and  3.  Then  Ki{h2{i))  =  an,o  +  an,2,  «2(A2(0)  =  l  +  2a„,i+a2  ^  +  20^  2  +  15a^^  + 

6a„,3  +  6a„,ia„,3,  «3(^2(0)  =  ^^n,2  +  12a„,ia„,2  +  72a„,2an,3  +  o(n''),  /C4(h2(0)  = 
24a„.3+48a„,ia„,3+48a^_2+432a^3+o(n~*).  If  we  need  |/c,(A:„)-/c,(/i2(<f))l  = 
for  each  t,  1  <  i  <  4,  the  unique  solution  of  a,’s  in  /i2(0^  order  0(n~^)  is 

such  that  a„,o  =  k.i  -  ln,z/Q  +  o(n"'),  a„,i  =  /„_2/2  -  +  (5/36)/^  3  +  o(n"'). 

an, 2  =  ^71,3/6  +  o(n"'),  and  a„,3  =  ln,4/'2i  -  ll^s/lS  +  o(n“^).  □ 

Proof  of  Proposition  3.1;  Results  (1),  (2),  and  (4)  follow  from  the  definition  of 
{F,  :  i  >  !}•  For  (3),  we  have 

M/3i)  +  iM/?2)  +  6/i(/33)  +  --- 

<k{ai)  +  [^(02)  +  •  •  •  +  /i(a6+i)]  4-  [^(06+2)  +  •  •  •  +  ^(q^26+i)]  +  ■  •  • 

t=i 

=  [/l(Oi)  +  •  •  •  +  /i(0(,)]  +  [/l(Qi+i)  + - 1"  ^(026)]  +  •  •  • 

<bh{0i)  +  hh[^2)  4-  hh[^3,)  +  •  •  • , 

where  the  first  inequality  comes  from  /?,  =  Q(i-i)bAi  <  Q(,_i)(,+i_j,  for  all  j  >  0 
and  the  monotonicity  of  h.  Similajly,  the  l^lst  inequality  is  a  consequence  of  /3,  > 
for  all  ;  >  1  and  the  monotonicity  of  h.  The  result  then  follows.  Q 

Two  more  prositions  remain  to  be  proved.  These  v,ill  require  a  number  of 
lemmas  which  are  given  below. 
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Lemma  A.  1.  Suppose  Y  €  (T(Xi,...,Xk)  and  Z  £  a(Xk+„,Xii+n+i,-- 

(1)  If\Y\  is  bounded  by  C  and  \Z\  is  bounded  by  D,  then 

|cum(y,  Z)|  <  4C£>a„. 

(2)  IfE[Y*]  <  C  and  E[Z*]  <  D,  then 

lcum(r,Z)l<4(l  +  C  +  I>)ay^ 

(3)  If  E[Y^]  <  C  and  E[Z*]  <  D,  then  for  positive  numbers  M  and  N, 

|cum(y,  Z)\  <  4.UA  a„  +  WNAr'^  +  4CMN-^  + 

Proof:  (1)  See  Lemma  2  of  Billingsley  [5],  p.  317. 

(2)  See  Lemma  3  of  Billingsley  (5],  p.  317. 

(3)  Let  A  be  a  set  and  define  the  indicator  function  1^  take  value  1  on  A  and  0  else¬ 

where.  Let  Yq  =  YI^y\<n,  Yi  =  F/|k|>n,  Zq  =  YI\z\<m,  and  Zi  =  ZI\z\>m.  Then 
\E[YoZo]-E[Yo]E[Zo]\<4MNa^,  \E[YoZ^]-E[Yo]E[Zi]\  <ADM-^N,  \E[Y,Zo]- 
E[Yi]E[ZQ]\  <  ACMN-^,  and  \E[YiZi]  -  E[Yi]E[Zi]\  <  The 

result  follows  from  lcum(y,  Z)\  =  |£:[yZ]  -  E[Y]E[Z]\  <  \E[YoZo]  -  E[Yo]E[Zo\\  + 
\E[YoZ,]  -  E[Yo]E[Z,]\ -h  \E[Y,Zo]  -  E[Yi]E[Zo]\  +  \E[Y,Z^]  -  E[Y,]E[Z,]\.  □ 

Lemma  A. 2.  Suppose  W  €  <7(Ai,...,A';t),  Y  £  a(Xjfc+„, . . . ,  A";),  and  Z  £ 
5 . . .).  Let  Qfn,n  —  min Then 

(1)  If\W\  is  bounded  by  B,  |y|  is  bounded  by  C,  and  \Z\  is  bounded  by  D,  then 

IcumflL',  y,  Z)|  <  SBC  Dam, n- 

(2)  HE[\W^\]  <  B,  E[\Y‘^\]  <  C,  and  E[\Z^\]  <  D,then 
|cum(iyy,Z)|  <  SLMNam,n 
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+  &MNBL-^  +  SLNCM-^  +  8LMDN~^ 

+  +  2MB^I'^D^I'^L-'^N-'^ 

+  +  B^/3(Ol/3pl/3^-l/3^^-l/3^-l/3 

(3)  UE[\W*\]  <  B,  ^[ly^l]  <  C,  and  E[\Z^\]  <  D,then 

\cnn^W,Y,Z)\  <  caifi, 

where  c  is  a  constant. 

Proof:  (1)  Notice  that  |cum(H';  KZ)|  =  -  £’(ir)][r  -  E(y)}[Z  -  E(Z)]I  < 

2B!EIV  -  E{Y)][Z  -  E[Z)]\  <  SBCDon.  Similarly,  lcum(Vr,  V', Z)1  <  8PCZ)q„ 
and  |cum(l'V,  K,  Z)|  <  SBCDam+n- 

(2)  Similar  to  the  proof  of  Lemma  A. 1(3). 

(3)  The  desired  result  can  be  obtained  by  using  (2)  and  let  Z-  =  A/  =  N  ~ 

□ 


Lemma  A. 3.  Suppose  <  oo- 

(1)  If  EYq  <  Cb~^  for  some  constant  C,  then  Yl'^=-oo  lCov(y'o, Ki)|  =  0{b~^). 

(2)  If  EYq  <  Db~*  for  some  constant  D,  then  Yl'n=-oo  ~  0{b~^). 

(3)  If  EYq  <  Db~*  for  some  constant  D,  then  X^^_oc  |Cov(y^^,  F„^)|  =  0(b~^). 

Proof:  (1)  From  Lemma  A.l,  |cum(yoi^t)|  £  +  4C6‘“'.V~^  +  \Cb~'^N~^  + 

Cb~^N~^.  Let  N  =  the  above  ineaqality  can  be  reduced  to  )cum(}o,  i’,)| 

<  cb-^3]'^  so  that  YT=-o^  |cum(ro,V;)|  <  cb'^  T.Z-oo^i^^  ^  = 

(2)  From  Liapounov  inequality,  follows  that  E'iQ  < 
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Cb  ^  where  C  =  From  Lemma  A.  1,  |cum(yo,  <  iMN^i+4:Db~‘^M~^N-\- 
ACb-\\fN-^  +  C^I^D^I'^b-^M-^N-\  Let  .V'  =  and  M  =  b~^3;^'\  the 

above  ineaqality  can  be  reduced  to  |^um(lo)  so  that 

E”-„|c''m(u,v;=)i  <  i  cb-^i-^YZ  But 

there  is  no  0{b~^^^)  terms  on  the  left  hand  side,  it  follows  that  J3^-oo  |Cov(yo,  Y3)\ 
is  at  most  0(6”^). 

(3)  From  Lemma  A.l,  |cum(y 0^,  y’'^)|  <  4A'^i9,  +  9Db~'*N~^.  Let  A'  = 

then  |cum(yo^,  }''■)!  <  cb~^  Yl^_^  3^^ ■  The  result  now  follows.  □ 

Lemma  A. 4.  Assume  lCm.n=-oo 

(1)  I{  EYq  <  Cb~^  for  some  conAants  C, .  then  X^m,n__oc  |cum(yo,  Ij,)!  = 
0(6-^). 

(2)  If  r^  'f  <  Db~*  for  some  coastaat  D,  then  ]C^,„=_oo  |cum(yc).y'',r.  VCf)]  = 
0(6-2). 

f3j  If  EYq  <  Db~^  for  some  constant  D,  then  X]m,n=-oo  = 

0(6-3). 

(4)  If  EYq  <  Db~*  for  some  constants  D.  then  X]m,n=-oo  )l  = 

0(6-3). 

Proof:  (1)  Let  L  =  b~^^^33i^n‘'.  M  =  b~^^^ 33^.1*  1  and  A'  =  b~^^'^3^,n^  and  then 

from  Lemma  .A. 2.  ]cum(io,  <  cb~^3mX^  where  c  is  a  constant,  so  that  the 

result  follows. 

(2)  Sirmlar  to  (1)  e.xcept  choosing  L  =  b~^^^ 3rn.l^ ,  -^f  =  b~^^^3m]i^.  ?nd  A’  = 
6-’ 


79 


(3)  Similar  to  (1)  except  choosing  L  =  b  M  =  b  and  N  = 

V  ,n  • 

(4)  Simileir  to  (1)  except  choosing  L  =  b~^0m]n*,  M  =  b~^l^m]n^,  and  N  =  b~^j3m]n*. 


Remark  A. 5.  Basically,  in  Lemmas  A. I  through  A. 4,  we  have  demonstrated 
that  an  infinite  sum  of  cumulants  such  as 

can  be  calculated  by  the  following  approach.  Ue  first  show  that  each  summand 
is  bounded  by,  say,  cb~'h{l3,ni, . . .  ,nk),  where  c  is  a  constant,  b  is  the  batch  size, 
and  h(-)  is  a  real  value  function.  However,  by  suitable  assumptions  of  bounded 
moments  and  asymptotic  properties  of  mixing  constants,  the  infinite  sum 


'^h{d,ni,...,nk) 

is  finite  so  that  ^  |cum[y 0° ,  , . . . ,  is  of  order  0{b~'). 

Corollary  A. 6  (Cumulants  of  A„).  Assume  that  {X,  :  i  >  1}  is  a  discrete 
time,  strictly  stationary  stochastic  process  with  zero  mean. 

(1)  If  £’lXil'‘  <  00,  and  the  sequence  is  mixing  with  q„  =  for  some  e  >  0, 

then  /ci(X^)  =  C>(m-^). 

(2)  If  £1X1!'^  <  00,  and  the  sequence  is  mixing  with  for  some  e  >  0. 

then  KiiXl)  =  0{m-^). 

{■3J  If  EjX I  <  oc,  and  the  sequence  is  mixing  with  q„  =  for  some 

c  >  0.  then  K3(.\'1,)  =  0(m~^). 
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(4)  If  E\Xi\^  <  oo,  cLnd  the  sequence  is  mixing  with  a„  =  for  some 

£  >  0,  then  K4(X^)  =  0{m~*). 

Proof:  All  four  results  are  special  cases  of  Lemma  2.6.  Q 

Corollary  A. 7  (Mixed  cumulants  of  and  A^).  Assume  that  {A,  :  t  >  1} 
is  a  discrete  time,  strictly  stationary  stochastic  process  with  zero  mean. 

(1)  If  E\Xi  I*  <  oo,  and  the  sequence  is  mixing  with  On  =  for  some  e  >  0, 

then 

(2)  If  E\Xi  <  DC,  and  the  sequence  is  mixing  with  a„  =  for  some  e  >  0, 

then  K2,i(^r.Xj  = 

(3)  If  £|A'’ip®  <  oo,  and  the  sequence  is  mixing  with  q„  =  for  some 

€  >  0,  then  Ki,2(Xr„,Xl)  =  and  K3,i(A,n,T^)  -  0(m-^). 

(4)  If  E\Xi\‘^  <  00,  and  the  sequence  is  mixing  with  a„  =  for  some 

e  >  0,  then  K2,2(A,n,  A^)  =  0(m"^)  and  «:4,i(Am,  A^)  =  0{m~*). 

(5)  If£|Ail^‘'  the  sequence  is  mixing  with  q„  =  0(n~^‘*~‘)  forsomee  >  0, 

then  ki.3(A„,A^)  =  0(m-'‘),  K3,2{A,„,A^)  =  0(771-“),  and  Ks,i{Xm,xlj  = 
0(m-5). 

(6)  if  £lAi|^  <  oo,  and  the  sequence  is  mixing  with  a„  =  0(n~^^~‘)  forsomee  >  0, 

then  «2,3(Am,A^)  =  0(m"“),  /C4,2(Am,A^)  -  0(m“®),  and  K6,i(A,„,A^)  = 
0(m-«). 

(7)  If  ElXil^"^  <  oc,  and  the  sequence  is  mixing  with  q„  =  0(n~^^~‘)  for  some  e  >  0. 

then  =  0(m-*),  «5.2(A,n,A^)  =  0(771-"),  and  «7.i(A,„,  A^)  = 

0(771-’;. 
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Proof;  These  results  can  be  obtained  from  Lemma  2.6.  [~] 

Lemma  A. 8  (Cumulants  of  (1/n)  Assume  that  {-Y,  :  i  >  1}  is  a 

discrete  time,  strictly  stationary  stochastic  process  with  zero  mean. 

(Ij  If  £^1-Yi|^  <  oo,  and  the  sequence  is  mixing  with  an  =  0(n~^),  then 

(2)  If  E\Xi^'^  <  oo,  and  the  sequence  is  mixing  with  a„  =  0(n“®),  then 

(3)  If  £’|A^ip°  <  oo,  and  the  sequence  is  mixing  with  a„  =  0{n~^^),  then 
«3((l/n)Er=i>^^)  =  0(n-W")- 

(4)  If  E\Xi^^  <  oo,  and  the  sequence  is  mixing  with  an  =  0(n~^'^),  then 
«4(l/n)Er=i>^^)  =  C(n-^)0(6-‘‘). 

Proof:  (1)  «x((l/n)  C=i  =  «2(Kx)  =  0(b-^). 

(2)  From  Lemma  2.8,  the  most  significant  term  of  K2((l/n)  J2?=i 

( V”)  which  is  (l/n)0(6“^)  form  Lemma  A. 2. 

(3)  From  Lemma  2.8,  the  most  significant  term  of  «3((l/n)  is 

(1/n^)  cum(y'o^,  y^^),  which  is  (l/n*)0(6“^)  from  Lemma  A. 2. 

(4)  This  can  be  proved  in  a  similar  fcishion  as  (2)  and  (3).  Q 

Lemma  A. 9  (Mixed  cumulants  of  X^  and  ' yj^)-  Assume  that 
{Xi  :  i  >  1}  is  a  discret  -  time,  strictly  stationary  stochastic  process  with  zero 
mean. 

(1)  If  E|Ai|®  <  oo,  and  the  sequence  is  mixing  with  an  =  0(n“'),  then 
^uiXmAW  LI:  y^)  =  0(n-^)0(b-^). 

(2)  If  <  oc,  and  the  sequence  is  mixing  with  q„  =  0(n~^),  then 
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«w(7„,(l/")ELi  y?)  =  0(n-^)0(b-^). 

(3)  If  E\Xi\^^  <  oo,  and  the  sequence  is  mixing  with  On  =  0(n~^^),  then 

ELi  y?)  =  0(n-')0(b-%  aod  K3.i(X„,  (1/n)  E",,  Y})  = 

0(n-S)0(6-3). 

(4)  If  E\Xif°  <  oo,  and  the  seouence  is  mixing  with  a„  =  0(n~^^),  then 

«2.2(X,„,  (1/n)  E?=i  )  =  0(n-3)0(6-3)  and  (1/n)  E"=i  = 

0(n-^}0(b-^). 

(5)  If  E\Xi\‘^‘^  <  oo,  and  the  sequence  is  mixing  with  a„  =  0(n"‘®),  tien 

/ci,3(X,„,(l/n)ELi  =  0(n-^)0(6-^),  /C3,2(X^,(l/n)ELi  i?)  = 

0(n-‘)0(6-^),  and  /C5,i(^m,  (1/n)  Er=i  ^")  =  0(n-^)0(6-*). 

(6)  If  <  oo,  and  the  sequence  is  mixing  with  a„  =  0{n~^^),  then 

K2,3(X,n,  (1/n)  Er=i  y,^)  =  0(n-^)0(b-%  (1/n)  Eh  = 

0(n-^)0(b-%  /C3,3(X,„,  (1/n)  E?=i  >^")  =  0(n-5)0(6-5),  and 
«6.i(X.,  (1/n)  E:=i  V;2)  =  0(n-«)0(6-^). 

(7)  If  ElXif^  <  oo,  and  the  sequence  is  mixing  with  o„  =  0{n~^^),  then 

«s.2(X^,  (1/n)  ELi  >;^)  =  0(n-^)0(b-^)  and  xr.j(X^,  (1/n)  Eh  Yi")  = 
0(n-")0(6-"). 

Proof:  (1)  From  Lemma  2.8,  the  most  significant  term  of  Ki^i(Xm,  (1/n)  Eh  ^i‘) 
is  (1/n)  E^-oo  ^^^(^0’ ^^)’  which  is  {l/n)0{b~'^)  form  Lemma  A. 3. 

(2)  From  Lemma  2.8,  the  most  significant  term  of  «:2,i(Xm,  (1/n)  Eh  Yi^)  i® 

(1/n^)  E~=_oo  which  is  (l/n’^)0(b~'^)  form  Lemma  A. 4. 

(3)  From  Lenuna  2.8,  the  most  significant  term  of  Ki,2(A,„,  (1/n)  Eh  Y,^)  is 
(1/n^)  Er^=_co  ^'^^(^0’ which  is  (l/n^'iO(6~^)  form  Lemma  A. 4, 

(4) -(7)  Can  be  shown  by  using  the  method  specified  in  Remark  A. 5.  E 


83 


Lemma  A. 10  (Mixed  cumulants  of  and  (V”)  2^"=!  Assume  that 

{X,-  :  i  >  1}  is  a  discrete  time,  strictly  stationary  stochastic  process  with  zero 
mean. 

(1)  If  <  oo,  and  the  sequence  is  mixing  with  =  0(n“®),  then 

=  0(n-=)0(6-2). 

(2)  If  £1X1!^°  <  00,  and  the  sequence  is  mixing  with  a„  =  0(n"'^),  then 

K^jxl  (1/n)  Zli  y?)  =  0(n-^)0(b-^)  and  (1/n)  ZL:  = 

0(n-3)0(6-3)  +  0(n-2)0(6-4). 


Proof;  (1)  Notice  that  (1/n)  y;^)  =  K2.i(A"m,  (1/n)  F,^).  The 

order  can  be  obtained  from  Lemma  A. 9(2). 


(2)  It  can  be  shown  that  (1/n)  V'“)  + 

4a3(X„)/t,,,(X„, (1/n)  r^)  +  4«,(X„)aj,,(X„,  (1/n)  E",,  V/).  The  desired 

order  follows  from  Lemma  A. 9.  The  second  result  can  be  proved  by  noting 


(1/")  ELi  >;”)  =  (1/n)  Er=.  r')+24,(X„,  (1/n)  ZL  and 


by  another  application  of  Lemma  A.9. 


□ 


_ 

Lemma  A. 11  (Mixed  cumulants  of  A'n,»  and  (1/n)  Y,^).  Assume 

that  {X,  :  i  >  1}  is  a  discrete  time,  strictly  stationary  stochastic  process  with  zero 
mean. 

(1)  If  L’lXip®  <  00,  and  the  sequence  is  mixing  with  a„  =  0(n~^^),  then 

(1/n)  ELi  i?)  =  0(n-^)0ib-^). 

(2)  if  £|Xip°  <  00,  and  the  sequence  is  mixing  with  a„  =  0(n~^^),  then 
nr...i(X..,Xl,(l/n)2/.,  1'4)  =  0(n-^)0(b-’). 

(3)  If  E\X\\^^  <  DC,  and  the  sequence  is  mixing  with  q„  =  0(n“^®),  then 


84 


0(n-3)0(6-‘‘),  and  /C3,i,i(^->^L  (1/n)  ^^=1  =  0{n-*)0{b-*). 

(4)  li <  00,  and  the  sequence  is  mixing  with  q„  =  0(n~^^),  then 

K2.2,i(Xrn,xlil/n)Z:^,y;^)  =  0{n-<)0ib-^),  K,,^,,(X^,xUlln)Z7=r  = 

0in-*)0{b-^),  and  (1/n)  ELi  V;")  =  0{rr^)0{b-^). 

(5)  If  <  oo,  and  the  sequence  is  mixing  with  q„  =  0(n“*^),  then 

K,,2AXm,xl^  (1/n)  ELi  i?)  =  Oin-^)0{b-^)  and  K,,,/,X^,xlAl/n)  Er=i 

=  0(n-®)0(6-"). 

Proof:  For  each  of  the  desired  results,  we  first  derive  an  identity  and  then  derive 
the  order  from  Lemma  A. 9.  Those  identities  are  as  follows. 

(1)  xL,  (1/n)  E"=i  =  «3.,(X,„,  (1/n)  ELi  + 
2«2(X,„)«i,i(^m,(l/n)E,’Li  y^)- 

(2)  K2.u(X,„,X1,  (1/n)  Er=i  =  «4.i(X„.,  (1/n)  E"=i  Y^  + 

2Ks(Xm)KiAXm,  (1/n)  Er=l  Y,^)  +  i>^2(Xm)K2A(Xm,  (1/n)  ELl  Y,^)' 

(3)  «i.2,i(X,„,  xl,  (1/n)  Er=i  1".^)  =  «5.i(X,„,  (1/n)  ZU  >?)  + 

4/t4(X,„)«i.i(X„„  (1/n)  Er=i  1?)  +  8«2(X,„)/C3.i(X„.,  (1/n)  Er=i  Y,^)  + 
8«3(X,„)k2.i(X„.,  (1/n)  ELi  V;2)  +  8k1(X,„)«,.i(X„.,  (1/n)  Er=i  Y^, 
Ka.i,2(X^,xl,  (1/n)  E?=x  Y^)  =  «3.2(X,„,  (1/n)  Er=i  Y,^) 

+  4«2.i(X„.,  (1/n)  ELi  i?)Ki,i(X,„,  (1/n)  ELi  + 

2«2(X„.)«i,2(X^,(l/n)  E:=i  V?);  and 

K3.i,i(X,„,xL,  (1/n)  E:=i  Yn  =-  'Cs.ilX^^d/n)  E^li  Y,^)  + 

2«4(X,„)«i,i(X„.  (1/n)  ELi  1?)  +  6«2(X,n)K3.i(X„.,  (1/n)  E?=i  Y,^)  + 

6«3(X,„)«2.i(Xm.(l/n)E?=i>';')- 
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(4)  «2.w(X„,TL,(l/n)£;L,y'4)  =  -C6..{X„,(l/n)Er.,i^') 

+  12«j(X„)«,,,(X„,(l/")Er=>  y?)  +  4/C5(X„)/c,,i(X„,(l/n)  + 

16/C3(X„)«3,.(X„,  (1/n)  X;,'L,  5^=)  +  12<C4(X„)«„(X„,  (l/n)  V'^)  + 

24«=(X„)k3,.(X„,  (l/n)  y;=)  +  24<C3(X„),cj(X„)K3,,(X„,(l/n)  j;”,, 

n2,i.!(X„,Xl,  (l/n)  i:”,,  y;“)  =  n,,3(X„,  (l/n)  + 

4«3.i(X„,(l/n)  J:",,  (l/n)  el,  Y?)  + 

2n3(X„)/c.,,(X„,  (l/n)  4/’)  +  4/C3(X„)n3.3(X„,  (l/n)  V;:')  + 

44,(X„,(l/n)EL,4^');  ^n<i 

n,,..i(X„,Xi,(i/n)i;/,.  Y^)  =  /<6.,{X„,(i/n)£"„,  y;^)  + 

8/C2(X„)n„(X„,  (l/n)  el,  K')  +  2«3(X„)«,,,(X„,(l/n)  EL,  V)  + 

12n3(X„)/C3.,(X™,  (l/n)  EL,  4^*)  +  8»<(^-)nM(X„,  (l/n)  el,  4',")- 

(5)  /C3.2.,(X„,XL  (l/n)  el,  4^^)  =  «T,,(X„,  (l/n)  EL,  4(4)  + 

16n2(X„)n5.,(X„.  (l/n)  EL,  4(4)  +  4n6(X„)n,,,(X„,(l/n)  EL,  Y?)  + 
28n3(X„)n4,,(X„,  (l/n)  EL,  Y?)  +  16ni(X„)n2,,(X„,  (l/n)  EL,  Y?)  + 
28/t,(X„)n3,,(X„,  (l/n)  EL,  4(4)  +  484(X„)n3.,(^m,  (l/n)  EL,  Y^)  + 
244(X„)/<,,,(X„,(l/n)EL,4(4)+32n4(X„)n,(X„)n,,,(X„,(l/n)EL,4’4)  + 
96/C3(X„)k2(X„)/<2,i(J1„,  (l/n)  el,  4(4);  nnd 
n3,,,2{X„,Xl,(l/n)EL.  4(4)  =  «3.2(X„,(l/n)  EL,  4(4)  + 

6»2(X„)«3.2(X„,  (l/n)  el,  '"4)  +  6K3(X„)n2.2(X„,(l/n)  el,  4(4)  + 
4n,..(X„,(l/n)  el,  y(4)n,,,(X„,(l/n)  EL.  4(')  + 

12n3,,(^,».(l/")  EL.  i(4)n2,.(^,n,(l/n)  EL,  4(4)  + 

2«4(X„)«,,2(X„,  (l/n)  el.  4(4)  + 

3/<J(X„)/c,,2(X„,(l/n)EL,4;4).  □ 
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Lemma  A. 12  (Cumulants  of  A„,i).  Assume  that  {X,  :  i  >  1}  is  a  discrete 
time,  strictly  stationary  stochastic  process  with  zero  mean. 

(1)  If  E\Xi  1^  <  oo,  and  the  sequence  is  mixing  with  a„  =  then  /Ci(A„,6)  = 

0{b-^)  +  0{n-^) 

(2)  If  E\Xi\^^  <  oc,  and  the  sequence  is  mixing  with  a„  =  0{n~^),  then  K2(A„,i)  = 

(3)  If  E\Xi\'^'^  <  oc,  and  the  sequence  is  mixing  with  then  «3(A„,fc)  = 

0(n-^) 

Proof:  RecaU  that  A„,i  =  V^j.l[^llb)  -  1  and  =  (1/n)  ¥/  -  xi,  then 

(1)  =  E[V„ilWilb)  -  1|  =  (6/oi)£((l/n)  2^=.  -  D  "  WOE'Til 

=  (blcl)E(Yf  -  1)  -  (blai)E-Ji  =  0(6-)  +  O(n-). 

(2)  =  (6/:7i)V4(I/n)  el.  1?  -  =  (6/,7=  )^[0(n-)0(6-)  + 

0(n“^)0(6"^)  +  0{n~'^)0{b~"^)]  =  0{n~^),  where  the  orders  are  obtained  from 
Corollary  A. 6,  Lemma  A. 8,  and  Lemma  A.  10. 

(3)  «3(An.i)  =  =  {b/crin0{n-^)0{b-^)]  =  0{n~% 

where  the  orders  are  obtained  from  Corollary  A. 6,  Lemma  A. 8,  and  Lemma  A. 10. 

□ 

Lemma  A. 13  (Mixed  cumulants  of  and  A„,i,).  Assume  that  {A,  ;  i  >  1} 
is  a  discrete  time,  strictly  stationary  stochastic  process  with  zero  mean. 

(1)  If  E\X  if  <  oo,  and  the  sequence  is  mixing  with  Qn  =  0(n“^),  then 
/cj.i(A„.i,A,„)  =  0(n-‘)0(6-^). 

(2)  If  Alp®  <  30,  and  the  sequence  is  mixing  with  q„  =  then 
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(3)  If  <  CO,  and  the  sequence  is  mixing  with  a„  =  0(n  ®),  then 

K^,,{A,,,,Xrr.)  =  0in-^)0(b-^). 

(4)  U  E\Xi'^*  <  oo,  and  the  sequence  is  mixing  with  a„  =  0(n~^^),  then 

(5)  If  <  oc,  and  the  sequence  is  mixing  with  q„  =  then 

«2.2(A„,i,X,„)  =  0(n-3)0(6-^). 

(6)  If  EjXij^^  <  oo,  and  the  sequence  is  mixing  with  a„  =  then 

«i,3(A„,6,X,„)  =  0(n-3)0(6-2). 

(7)  If  E\Xi\‘^^  <  oo,  and  the  sequence  is  mixing  with  a„  =  C>(n"^'),  then 
«3,2(An,6,r^)  =  6>(n-^)0(A-l). 

(8)  If  E\Xi\‘^*  <  oc,  and  the  sequence  is  mixing  with  a„  =  0(n“^®),  then 

«2.3(A„,6,X,„)  =  0(n-'‘)C»(6-2). 

(9)  If  E\Xi\^^  <  00,  and  the  sequence  is  mixing  with  a„  =  0(n~^^),  then 
.KU^„,,,Xm)  ^  0{n-*)0{b-^). 

(10)  /f  E|A'ip^  <  oo,  and  the  sequence  is  mixing  with  q„  =  C>(n~'®),  then 
/C3,3(A„,i,A,„)  =  0(72-5)0(6-2). 

(11)  If  £'|Aip®  <  oo,  and  the  sequence  is  mixing  with  q„  =  0(77“^^),  then 
«2,4(A„,fc,X,„)  =  0(71-5)0(6-3). 

Proof:  RecaU  that  A„,,  =  K,.6/(<rJ./6)  -  I  and  K,.,  =  (1/n)  then 

(1)  =  (V<ti)[at,i((l/n)  EL,  +  k.,.(X1.X„)1  = 

0(71-')0(6-^).  where  the  orders  are  obtained  from  Corollary  A. 7  and  Lemma  A. 9. 

(2)  a2..(A„,..X„)  =  (6/<^L)"*,..,iI(1/")EL.7;=  -Xl,(l/n) 

=  0(72-^)0(6-'),  where  the  orders  are  obtained  from  Corollary  A. 7,  Lemma  A. 9 
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and  Lemma  A.ll. 


(3)  -  xl,X„,X„]  =  0(n-')0(i-), 
where  the  orders  are  obtained  from  Corollary  A. 7  and  Lemma  A. 9. 

(4) 

=  (6/cr^)^0(n“^)C)(6“‘‘)  =  0{n~^b~^),  where  the  orders  are  obtained  from  Corol¬ 
lary  A. 7,  Lemma  A. 9,  and  Lemma  A.ll. 

(5)  K2.2(A„.fc,X„)  =  (6/a^J^K2.2l(l/n)£r=i  -XL,X„] 

~  (6/o‘^)^0(n“^)0(6“^)  =  0in~^)0{b~^),  where  the  orders  are  obtained  from 
Corollary  A. 7,  Lemma  .A.. 9.  and  Lemma  .A.ll. 

(6)  ,C.,3(A,,s,A„)  =  (V£'i)^-l.3l(l/")EL, 

—  =  0(n~^)0(b~^),  where  the  orders  are  obtained  from 

Corollary  A. 7,  Lemma  .A. 9,  and  Lemma  A.ll.- 

(7)  «3,2(a„..,X„)  =  (4M)’'C3,3i(i/-.)  r",,  y? 

=  {bf (T^)^0{n~‘*)0{b~*)  =  0(n~‘*)0(b~^),  where  the  orders  are  obtained  from 
Corollary  A. 7,  Lemma  A. 9,  and  Lemma  A.ll. 

(8)  /<3,3(A„,»,X„)  =  (WO‘«r3Ul/n)Z:.,  V;'  -XL,X„J 

=  {bj a‘^)‘^0{n~'*)0{b~^)  =  0(n~‘*)0(b~^),  where  the  orders  are  obtained  from 
Corollary  A. 7,  Lemma  A. 9,  and  Lemma  A.ll. 

(9)  Ki,4(An,b,X„)  =  ib/(Tl)Ki,4l(l/n)  J2?=l 

=  (b/<T^)^0(n~‘*)0(b~‘*)  =  0(n“‘‘)0(6“^),  where  the  orders  are  obtained  from 
Corollary  A. 7,  Lemma  A. 9,  and  Lemma  A.ll. 

(10)  K3,3(A„,6,X„)  =  (o/a^)^K3.3[(l/n))C.=i 

=  [b j al^)^0{n~^)0{b~^)  =  0(n~^)C>(6"^),  where  the  orders  are  obtained  from 
Corollary  A. 7,  Lemma  .A. 9,  and  Lemma  A.ll. 

(11)  K2,4(A„,6,X,„)  =  (6/c’m)^K2,4((l/«)  Er=i 


89 


=  (6/<7^)^0(n  ®)0(fe  =  0(n  ®)0(6  ^),  where  the  orders  are  obtained  from 

Corollary  A. 7,  Lemma  A. 9,  and  Lemma  A. 11.  Q 

Proof  of  Proposition  3.4:  Recall  that  {-Y,}  is  zero  mean.  From  equation  (3.10), 

<n,6  = 

{1  -  (l/2)A„.i  +  (3/8)Al,  -  (5/l6)Al^  +  •••}. 

(1)  =  (m^^ycr^)E[X^{l  -  (l/2)A„,t  +  (3/8)A2,,  -  (S/ie.^.A^,,  -  ••}]  = 

I (7rn)E{Xm)  —  (l/2)(m‘/^/(7m)£'[A'mA„,6]  +  •  •  ••  The  first  term  on  the  right 
hand  side  is  0,  the  second  term  is  and  has  the  most  significant 

term  — (l/2)(v/m6/<T^)/Ci,i(Am,  (1/n)  V'^),  and  the  remainder  terms  are  of 

(2)  K2{in,b)  =  (m/cr^)«2(Am)-(m/cT2  )/Ci,j(X,„,A„,fc)  +  (l/4)(m/c7^)AC2(An,6)+.... 
The  desired  result  then  follows  from  Lemmas  A. 12  and  A. 13. 

(3)  Note  that  =  {m^!^ IcI)k^(X,^)  -  (3/2)(m3/Vcr^)«2,,(X,„,  A^.j)  + 

(3/4)(m^/V<7m)«i.2('^m,  An,6)  “  ( l/8)(m3/^/<7^)«:3( A„  (,)  + - The  desired  result 

tehn  follows  from  Lemmas  A.  12  and  A.  13. 

(4)  Note  that  K3(tn,6)  =  (^Vo^r  '■"4(A,n)  ~2{m^la*^)K3,i[Xm,Aj,,h)  + 
(3/2)(mVcr;l,)«2,2(X,„,  A„.i)  -  (l/2)(mVcT:|.)«i.3(A,„,  A„.fc)  + 

(l/16)(m^/(T^)K4(A„,(,)  + _ The  desired  result  then  follows  from  Lemmas  A. 12 

and  A.13.  Q 
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